
# importing libraries
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import warnings

# Setting the warnings to be ignored
warnings.filterwarnings('ignore')


#market_data = pd.read_csv(r"/content/drive/MyDrive/FACT_MARKET_DEMAND_Ref.csv")
import pandas as pd
market_demand = pd.read_csv("FACT_MARKET_DEMAND-001.csv")


market_demand['CALORIC_SEGMENT'].fillna(market_demand['CALORIC_SEGMENT'].mode()[0], inplace=True)
print("\nMissing value imputation check:")
print(market_demand.isnull().sum().sum())

Missing value imputation check:
0


# filtering for Plum
plum = market_demand[(market_demand['ITEM'].str.contains('PLUM', case=False, regex=True))]


# Assuming 'DATE' column is in datetime format
plum['DATE'] = pd.to_datetime(plum['DATE'])

# Find the earliest date
earliest_date = plum['DATE'].min()

# Calculate the week numbers from the earliest date
plum['WEEK_NUMBER'] = (plum['DATE'] - earliest_date).dt.days // 7 + 1


# Aggregate data by week number
weekly_sales = plum.groupby('WEEK_NUMBER')['UNIT_SALES'].sum()

# Plot time series
plt.figure(figsize=(18, 6))
weekly_sales.plot(marker='o', color='blue', linestyle='-')
plt.title('Weekly Unit Sales')
plt.xlabel('Week Number')
plt.ylabel('Total Unit Sales')
plt.grid(True)
plt.show()


import pandas as pd
import matplotlib.pyplot as plt
from scipy.signal import argrelextrema

# Assuming 'plum' dataframe contains 'WEEK_NUMBER' and 'UNIT_SALES' columns

# Aggregate data by week number
weekly_sales = plum[plum['MANUFACTURER'] == 'SWIRE-CC'].groupby('WEEK_NUMBER')['UNIT_SALES'].sum()

# Find local maxima (peaks)
peaks_index = argrelextrema(weekly_sales.values, comparator=lambda x, y: x > y, order=5)[0]
peaks = weekly_sales.iloc[peaks_index]

# Plot time series with peaks
plt.figure(figsize=(17, 8))
weekly_sales.plot(marker='o', color='blue', linestyle='-', label='Total Unit Sales')
peaks.plot(marker='o', color='red', linestyle='', markersize=8, label='Peaks')

# Annotate peaks with week number
for week, sales in peaks.items():
    plt.annotate(f'Week {week}\n({sales} units)', xy=(week, sales), xytext=(-20, 20), textcoords='offset points',
                 arrowprops=dict(arrowstyle='->', color='black'))



plt.title('Weekly Unit Sales with Peaks')
plt.xlabel('Week Number')
plt.ylabel('Total Unit Sales')
plt.grid(True)
plt.legend()
plt.show()


plum['13_WEEK_INTERVAL'] = (plum['WEEK_NUMBER'] - 1) // 13 + 1


import matplotlib.pyplot as plt

# Assuming rolling_sales_df contains the DataFrame with 'SALES_13W_SUM' and '13_WEEK_INTERVAL' columns

# Extract data for plotting
sales_sum = rolling_sales_df['SALES_13W_SUM']
interval = rolling_sales_df['13_WEEK_INTERVAL']

# Plot the rolling sum
plt.figure(figsize=(18, 6))
plt.plot(interval, sales_sum, marker='o', linestyle='-')
plt.xlabel('13-Week Interval')
plt.ylabel('Sales Sum')
plt.title('Rolling Sum of Sales over 13-Week Intervals')
plt.xticks(rotation=45)  # Rotate x-axis labels for better readability
plt.grid(True)

# Find the index of the maximum value
max_index = sales_sum.idxmax()
max_value = sales_sum[max_index]
max_interval = interval[max_index]

# Annotate the peak with a red marker and the corresponding 13-week interval
plt.annotate(f'Max: {max_value} ({max_interval})',
             xy=(max_interval, max_value),
             xytext=(80, -20),
             textcoords='offset points',
             arrowprops=dict(facecolor='red', arrowstyle='->'),
             fontsize=10,
             color='red',
             ha='center')

plt.tight_layout()
plt.show()


import matplotlib.pyplot as plt

# Assuming rolling_sales_df contains the DataFrame with 'SALES_13W_SUM' and '13_WEEK_INTERVAL' columns

# Extract data for plotting
sales_sum = rolling_sales_df['SALES_13W_SUM']
interval = rolling_sales_df['13_WEEK_INTERVAL']

# Plot the rolling sum
plt.figure(figsize=(18, 6))
plt.plot(interval, sales_sum, marker='o', linestyle='-')
plt.xlabel('13-Week Interval')
plt.ylabel('Sales Sum')
plt.title('Rolling Sum of Sales over 13-Week Intervals')
plt.xticks(rotation=45)  # Rotate x-axis labels for better readability
plt.grid(True)

# Find the index of the maximum value
max_index = sales_sum.idxmax()
max_value = sales_sum[max_index]
max_interval = interval[max_index]

# Annotate the peak with a red marker and the corresponding 13-week interval
plt.annotate(f'Max: {max_value} ({max_interval})',
             xy=(max_interval, max_value),
             xytext=(80, -20),
             textcoords='offset points',
             arrowprops=dict(facecolor='red', arrowstyle='->'),
             fontsize=10,
             color='red',
             ha='center')

plt.tight_layout()
plt.show()


swire_plum = plum[plum['MANUFACTURER'] == 'SWIRE-CC']
swire_plum= swire_plum.sort_values('DATE')

swire_plum_daily = swire_plum.groupby('DATE')['UNIT_SALES'].sum().reset_index()
swire_plum_daily.head(2)


from prophet import Prophet
from prophet.plot import plot_plotly, plot_components_plotly
import pandas as pd


# Convert DATE column to datetime if not already in datetime format
swire_plum_daily['DATE'] = pd.to_datetime(swire_plum_daily['DATE'])

# Sort the data by date
swire_plum_daily = swire_plum_daily.sort_values(by='DATE')

# Calculate the index for partitioning the data
partition_index = int(len(swire_plum_daily) * 0.8)

# Split the data into train and test sets
train_data = swire_plum_daily.iloc[:partition_index]
test_data = swire_plum_daily.iloc[partition_index:]

# Rename columns as per Prophet's requirements
train_data = train_data.rename(columns={'DATE': 'ds', 'UNIT_SALES': 'y'})
test_data = test_data.rename(columns={'DATE': 'ds', 'UNIT_SALES': 'y'})

# Print the lengths of train and test dataframes
print("Length of Train Data:", len(train_data))
print("Length of Test Data:", len(test_data))

# Print the start and end dates of train and test dataframes
print("Start Date of Train Data:", train_data['ds'].min())
print("End Date of Train Data:", train_data['ds'].max())
print("Start Date of Test Data:", test_data['ds'].min())
print("End Date of Test Data:", test_data['ds'].max())

Length of Train Data: 119
Length of Test Data: 30
Start Date of Train Data: 2020-12-05 00:00:00
End Date of Train Data: 2023-03-11 00:00:00
Start Date of Test Data: 2023-03-18 00:00:00
End Date of Test Data: 2023-10-28 00:00:00


# Initialize and fit the Prophet model
p_model = Prophet()
p_model.fit(train_data)

# Create future dates for forecasting
future = p_model.make_future_dataframe(periods=34, freq='W')

# Make predictions
forecast = p_model.predict(future)

# Plot the forecast
plot_plotly(p_model, forecast)


plot_components_plotly(p_model, forecast)


pred = forecast[['ds','yhat']]


import matplotlib.pyplot as plt

# Plot test data
plt.figure(figsize=(12, 6))
plt.plot(train_data['ds'], train_data['y'], label='Actual Train Set', color='blue')

plt.plot(test_data['ds'], test_data['y'], label = 'Actual Test Set',color='lightblue')

# Plot forecast data
plt.plot(forecast['ds'], forecast['yhat'], label='Forecast', color='red')


# Add labels and title
plt.xlabel('Date')
plt.ylabel('Sales')
plt.title('Actual vs Forecast')
plt.legend()

# Rotate x-axis labels for better readability
plt.xticks(rotation=45)

# Show plot
plt.show()


# Function to calculate the scores

from sklearn.metrics import mean_absolute_error, mean_squared_error
from math import sqrt
import numpy as np
def get_scores(actual, predicted):
# Calculate errors
    mae = mean_absolute_error(actual, predicted)
    mse = mean_squared_error(actual, predicted)
    rmse = sqrt(mse)
    #percentage_diff = (np.abs(actual - predicted).sum()) /( np.abs(actual).mean())
    
    # Combine actual and predicted values into one DataFrame
    mape_df = pd.DataFrame({'actual': actual.values, 'predicted': predicted.values})

# Calculate absolute percentage error
    mape_df['absolute_percentage_error'] = np.abs((mape_df['actual'] - mape_df['predicted']) / mape_df['actual']) * 100

# Calculate MAPE
    mape = mape_df['absolute_percentage_error'].mean()
    
# Calculate MAPE
    #mape = percentage_diff.mean()
# Calculate "Accuracy" Percentage
    accuracy_percentage = 100 - mape
# Print metrics
#     print(f"MAE: {mae}")
#     print(f"MSE: {mse}")
#     print(f"RMSE: {rmse}")
#     print(f"MAPE: {mape}%")
#     print(f"Direct 'Accuracy' Percentage: {accuracy_percentage}%")
    return pd.Series(data={'mae':mae, 'mse':mse, 'rmse':rmse, 'mape':mape, 'direct_accuracy':accuracy_percentage}, index=['mae', 'mse', 'rmse','mape','direct_accuracy'])


from sklearn.metrics import mean_absolute_error, mean_squared_error
from math import sqrt
import numpy as np
import pandas as pd
from prophet import Prophet
import matplotlib.pyplot as plt

# Assuming mulberry_diet_sparkling_validation is your DataFrame
swire_plum_sparkling_validation = swire_plum_daily.copy()

# Convert Date to datetime variable
swire_plum_sparkling_validation['DATE'] = pd.to_datetime(swire_plum_sparkling_validation['DATE'])

# Group by date to ensure only one record for each date
swire_plum_sparkling_validation = swire_plum_sparkling_validation.groupby('DATE')['UNIT_SALES'].sum().reset_index()

# Sort the data by date
swire_plum_sparkling_validation = swire_plum_sparkling_validation.sort_values(by='DATE')

# Calculate the index for partitioning the data
partition_index = int(len(swire_plum_sparkling_validation) * 0.8)

# Split the data into train and test sets
train_data = swire_plum_sparkling_validation.iloc[:partition_index]
test_data = swire_plum_sparkling_validation.iloc[partition_index:]

# Renaming for Prophet requirements
train_data = train_data.rename(columns={'DATE': 'ds', 'UNIT_SALES': 'y'})
test_data = test_data.rename(columns={'DATE': 'ds', 'UNIT_SALES': 'y'})

p_model = Prophet()
p_model.fit(train_data)

# Create future dates for forecasting
future = p_model.make_future_dataframe(periods=0, freq='W')

# Make predictions
forecast = p_model.predict(future)

# Extracting the date and unit sales from forecast
pred = forecast[['ds','yhat']]

22:00:01 - cmdstanpy - INFO - Chain [1] start processing
22:00:01 - cmdstanpy - INFO - Chain [1] done processing


get_scores(test_data['y'], pred['yhat'])

mae                  162.088814
mse                43781.580074
rmse                 209.240484
mape                  13.020970
direct_accuracy       86.979030
dtype: float64


earliest_date = forecast['ds'].min()

# Step 2: Calculate week numbers
forecast['WEEK_NUMBER'] = ((forecast['ds'] - earliest_date).dt.days // 7) % 52
weekly_forecast_data = forecast.groupby('WEEK_NUMBER')['yhat'].sum().reset_index()


# Sort the DataFrame by 'WEEK_NUMBER' if it's not already sorted
weekly_forecast_data = weekly_forecast_data.sort_values(by='WEEK_NUMBER')

# Ensure that there are enough rows to compute the rolling sum
if len(weekly_forecast_data) < 13:
    print("Insufficient data to compute rolling sum.")
else:
    # Initialize an empty list to store the results
    rolling_sales_sum = []

    # Iterate over the DataFrame to calculate the rolling sum for each 13-week interval
    for i in range(len(weekly_forecast_data) - 12):
        # Calculate the sum of UNIT_SALES for the current 13-week interval
        sales_sum = weekly_forecast_data.loc[i:i+12, 'yhat'].sum()
        # Construct the 13-week interval string
        interval = f"{weekly_data.iloc[i]['WEEK_NUMBER']}-{weekly_data.iloc[i+12]['WEEK_NUMBER']}"
        # Append the result to the list
        rolling_sales_sum.append((sales_sum, interval))

    # Create a new DataFrame from the results
    rolling_sales_df = pd.DataFrame(rolling_sales_sum, columns=['SALES_13W_SUM', '13_WEEK_INTERVAL'])

    # Display the DataFrame with the rolling sum and 13-week intervals
# print(rolling_sales_df)


import matplotlib.pyplot as plt

# Assuming rolling_sales_df contains the DataFrame with 'SALES_13W_SUM' and '13_WEEK_INTERVAL' columns

# Extract data for plotting
sales_sum = rolling_sales_df['SALES_13W_SUM']
interval = rolling_sales_df['13_WEEK_INTERVAL']

# Plot the rolling sum
plt.figure(figsize=(18, 6))
plt.plot(interval, sales_sum, marker='o', linestyle='-')
plt.xlabel('13-Week Interval')
plt.ylabel('Sales Sum')
plt.title('Rolling Sum of Sales over 13-Week Intervals')
plt.xticks(rotation=45)  # Rotate x-axis labels for better readability
plt.grid(True)

# Find the index of the maximum value
max_index = sales_sum.idxmax()
max_value = round(sales_sum[max_index])
max_interval = interval[max_index]

# Annotate the peak with a red marker and the corresponding 13-week interval
plt.annotate(f'Max: {max_value} ({max_interval})',
             xy=(max_interval, max_value),
             xytext=(70, 20),
             textcoords='offset points',
             arrowprops=dict(facecolor='red', arrowstyle='->'),
             fontsize=10,
             color='red',
             ha='center')

plt.tight_layout()
plt.show()


plum_swire_diet = swire_plum[swire_plum['CALORIC_SEGMENT'] == 'DIET/LIGHT']
plum_swire_diet_weekly_sales = plum_swire_diet.groupby('DATE')['UNIT_SALES'].sum().reset_index()
plum_swire_diet_weekly_sales.head(2)


from prophet import Prophet
from prophet.plot import plot_plotly, plot_components_plotly
import pandas as pd


# Convert DATE column to datetime if not already in datetime format
plum_swire_diet_weekly_sales['DATE'] = pd.to_datetime(plum_swire_diet_weekly_sales['DATE'])

# Sort the data by date
plum_swire_diet_weekly_sales = plum_swire_diet_weekly_sales.sort_values(by='DATE')

plum_swire_diet_model_data = plum_swire_diet_weekly_sales.rename(columns={'DATE': 'ds', 'UNIT_SALES': 'y'})

# Initialize and fit the Prophet model
p_model = Prophet()
p_model.fit(plum_swire_diet_model_data)

# Create future dates for forecasting
future = p_model.make_future_dataframe(periods=52, freq='W')

# Make predictions
forecast = p_model.predict(future)

# Plot the forecast
plot_plotly(p_model, forecast)


earliest_date = forecast['ds'].min()

# Step 2: Calculate week numbers
forecast['WEEK_NUMBER'] = ((forecast['ds'] - earliest_date).dt.days // 7) % 52
weekly_forecast_data = forecast.groupby('WEEK_NUMBER')['yhat'].sum().reset_index()

# Step 4: Create a time series plot
plt.figure(figsize=(12, 4))
plt.plot(weekly_forecast_data['WEEK_NUMBER'], weekly_forecast_data['yhat'], marker='o', linestyle='-')
plt.xlabel('Week Number')
plt.ylabel('Unit Sales')
plt.title('Weekly Total Unit Sales')
plt.grid(True)
plt.show()


# Sort the DataFrame by 'WEEK_NUMBER' if it's not already sorted
weekly_forecast_data = weekly_forecast_data.sort_values(by='WEEK_NUMBER')

# Ensure that there are enough rows to compute the rolling sum
if len(weekly_forecast_data) < 13:
    print("Insufficient data to compute rolling sum.")
else:
    # Initialize an empty list to store the results
    rolling_sales_sum = []

    # Iterate over the DataFrame to calculate the rolling sum for each 13-week interval
    for i in range(len(weekly_forecast_data) - 12):
        # Calculate the sum of UNIT_SALES for the current 13-week interval
        sales_sum = weekly_forecast_data.loc[i:i+12, 'yhat'].sum()
        # Construct the 13-week interval string
        interval = f"{weekly_data.iloc[i]['WEEK_NUMBER']}-{weekly_data.iloc[i+12]['WEEK_NUMBER']}"
        # Append the result to the list
        rolling_sales_sum.append((sales_sum, interval))

    # Create a new DataFrame from the results
    rolling_sales_df = pd.DataFrame(rolling_sales_sum, columns=['SALES_13W_SUM', '13_WEEK_INTERVAL'])

    # Display the DataFrame with the rolling sum and 13-week intervals
    #print(rolling_sales_df)


import matplotlib.pyplot as plt

# Assuming rolling_sales_df contains the DataFrame with 'SALES_13W_SUM' and '13_WEEK_INTERVAL' columns

# Extract data for plotting
sales_sum = rolling_sales_df['SALES_13W_SUM']
interval = rolling_sales_df['13_WEEK_INTERVAL']

# Plot the rolling sum
plt.figure(figsize=(18, 6))
plt.plot(interval, sales_sum, marker='o', linestyle='-')
plt.xlabel('13-Week Interval')
plt.ylabel('Sales Sum')
plt.title('Rolling Sum of Sales over 13-Week Intervals')
plt.xticks(rotation=45)  # Rotate x-axis labels for better readability
plt.grid(True)

# Find the index of the maximum value
max_index = sales_sum.idxmax()
max_value = round(sales_sum[max_index])
max_interval = interval[max_index]

# Annotate the peak with a red marker and the corresponding 13-week interval
plt.annotate(f'Max: {max_value} ({max_interval})',
             xy=(max_interval, max_value),
             xytext=(70, 20),
             textcoords='offset points',
             arrowprops=dict(facecolor='red', arrowstyle='->'),
             fontsize=10,
             color='red',
             ha='center')

plt.tight_layout()
plt.show()


# week wise sales for best 13 weeks
best_13_weeks_q1 = forecast[(forecast['WEEK_NUMBER'] <= 12) & (forecast['ds'] > '2023-02-25')]
best_13_weeks_q1.head(2)


# Compute total sales for all 13 weeks
total_sales = best_13_weeks_q1['yhat'].sum()

# Compute lower and upper bounds for estimate
lower_bound = best_13_weeks_q1['yhat_lower'].sum()
upper_bound = best_13_weeks_q1['yhat_upper'].sum()

# Calculate confidence interval
confidence_interval = ((upper_bound - lower_bound) / total_sales) * 100

# Display the results
print("Total Sales for 13 weeks:", total_sales)
print("Lower Bound Estimate:", lower_bound)
print("Upper Bound Estimate:", upper_bound)
print("Confidence Interval (%):", confidence_interval)

Total Sales for 13 weeks: 9057.187762419222
Lower Bound Estimate: 6291.951059331268
Upper Bound Estimate: 11799.518942786954
Confidence Interval (%): 60.80880763351411


import matplotlib.pyplot as plt

# Extract week numbers and corresponding sales estimates
week_numbers = best_13_weeks_q1['WEEK_NUMBER']
sales_estimates = best_13_weeks_q1['yhat']
lower_bounds = best_13_weeks_q1['yhat_lower']
upper_bounds = best_13_weeks_q1['yhat_upper']

# Plot the time series graph with confidence interval
plt.figure(figsize=(12, 6))
plt.plot(week_numbers, sales_estimates, color='blue', label='Sales Estimate')
plt.fill_between(week_numbers, lower_bounds, upper_bounds, color='skyblue', alpha=0.4, label='Confidence Interval')
plt.xlabel('Week Number')
plt.ylabel('Total Sales')
plt.title('Week-wise Distribution of Total Sales with Confidence Interval')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()


plum_swire_diet['CATEGORY'].value_counts(normalize=True)

ENERGY             0.998592
SPARKLING WATER    0.001126
SSD                0.000282
Name: CATEGORY, dtype: float64


# Filter DataFrame for 'SSD' category
ssd_sales = market_demand[market_demand['CATEGORY'] == 'SSD']['UNIT_SALES'].sum()

# Filter DataFrame for 'ENERGY' category
energy_sales = market_demand[market_demand['CATEGORY'] == 'ENERGY']['UNIT_SALES'].sum()

print(" % of unit sales for SSD category compared to SSD + ENERGY category:", ssd_sales/(energy_sales+ssd_sales))
print(" % of unit sales for ENERGY category compared to SSD + ENERGY category:", energy_sales/(energy_sales+ssd_sales))

 % of unit sales for SSD category compared to SSD + ENERGY category: 0.6923682683459191
 % of unit sales for ENERGY category compared to SSD + ENERGY category: 0.30763173165408086


import pandas as pd

# Filter DataFrame for 'SSD' category
ssd_data = market_demand[market_demand['CATEGORY'] == 'SSD']['UNIT_SALES'].describe()

# Filter DataFrame for 'ENERGY' category
energy_data = market_demand[market_demand['CATEGORY'] == 'ENERGY']['UNIT_SALES'].describe()

# Concatenate the DataFrames horizontally
combined_data = pd.concat([ssd_data, energy_data], axis=1)
combined_data.columns = ['SSD', 'ENERGY']  # Rename columns for clarity


combined_data


plum_swire_diet['BRAND'].value_counts(normalize=True)

VENOMOUS BLAST         0.998592
CUPADA ARID            0.001126
DIET MITE PURE ZERO    0.000282
Name: BRAND, dtype: float64


# Filter DataFrame for 'SSD' category
vb_brand_sales = market_demand[market_demand['BRAND']=='VENOMOUS BLAST']['UNIT_SALES'].sum()

# Filter DataFrame for 'ENERGY' category
dietsmash_sales = market_demand[market_demand['BRAND']=='DIET SMASH']['UNIT_SALES'].sum()

print(" % of unit sales for DIET SMASH Brand compared to other two:", dietsmash_sales/(vb_brand_sales+dietsmash_sales))
print(" % of unit sales for VENOMOUS BLAST Brand compared to other two:", vb_brand_sales/(vb_brand_sales+dietsmash_sales))

 % of unit sales for DIET SMASH Brand compared to other two: 0.11897748583594121
 % of unit sales for VENOMOUS BLAST Brand compared to other two: 0.8810225141640589


# Filter DataFrame for 'SSD' category
vb_data =market_demand[market_demand['BRAND']=='VENOMOUS BLAST']['UNIT_SALES'].describe()

dietsmash_data= market_demand[market_demand['BRAND']=='DIET SMASH']['UNIT_SALES'].describe()

# Concatenate the DataFrames horizontally
combined_data = pd.concat([dietsmash_data, vb_data], axis=1)
combined_data.columns = ['DIET SMASH', 'VENOMOUS BLAST']  # Rename columns for clarity

combined_data


import matplotlib.pyplot as plt

# Filter the data for 'VENOMOUS BLAST' and 'DIET SMASH' brands
venomous_blast_data = market_demand[market_demand['BRAND'] == 'VENOMOUS BLAST']
diet_smash_data = market_demand[market_demand['BRAND'] == 'DIET SMASH']

# Create subplots
fig, axes = plt.subplots(1, 2, figsize=(12, 6))

# Plot for 'VENOMOUS BLAST' brand
axes[0].hist(venomous_blast_data['UNIT_SALES'], bins=30, color='skyblue', edgecolor='black')
axes[0].set_title('Unit Sales Distribution for VENOMOUS BLAST Brand')
axes[0].set_xlabel('Unit Sales')
axes[0].set_ylabel('Frequency')
axes[0].grid(True)

# Plot for 'DIET SMASH' brand
axes[1].hist(diet_smash_data['UNIT_SALES'], bins=30, color='lightgreen', edgecolor='black')
axes[1].set_title('Unit Sales Distribution for DIET SMASH Brand')
axes[1].set_xlabel('Unit Sales')
axes[1].set_ylabel('Frequency')
axes[1].grid(True)

# Adjust layout
plt.tight_layout()
plt.show()


import numpy as np
venomous_blast_data = market_demand[market_demand['BRAND'] == 'VENOMOUS BLAST']

#  percentile values
percentiles = np.arange(75, 96, 5)  # Percentiles from 75 to 95 in steps of 5
percentiles = np.append(percentiles, np.arange(96, 101, 1))  # Percentiles from 96 to 100 in steps of 1
percentile_values = np.percentile(venomous_blast_data['UNIT_SALES'], percentiles)

percentile_table = pd.DataFrame({'Percentile': percentiles, 'Value': percentile_values})
percentile_table


avocado = market_demand[(market_demand['ITEM'].str.contains('AVOCADO', case=False, regex=True))]
avocado_regular = avocado[(avocado['CALORIC_SEGMENT'].str.contains('REGULAR', case=False, regex=True, na=False))]
avocado_regular_ssd = avocado_regular[(avocado_regular['CATEGORY'].str.contains('SSD', case=False, regex=True))]
avocado_regular_ssd_swire = avocado_regular_ssd[(avocado_regular_ssd['MANUFACTURER'].str.contains('SWIRE', case=False, regex=True))]
avocado_regular_ssd_swire_small = avocado_regular_ssd_swire[(avocado_regular_ssd_swire['PACKAGE'].str.contains('SMALL', case=False, regex=True))]


df = avocado_regular_ssd_swire_small[['DATE','UNIT_SALES']].groupby(by="DATE", as_index=False).sum()
df['DATE'] = pd.to_datetime(df['DATE'])
df.set_index('DATE', inplace=True)
df = df.asfreq('W-SAT')


plt.plot(df)
plt.xlabel('Date')
plt.ylabel('UNIT_SALES')

Text(0, 0.5, 'UNIT_SALES')


# Train and Test periods
from datetime import datetime

s_date = '2023-02-25'
e_date = '2023-08-26'
df = df.loc[df.index <= e_date]
train_period = df.loc[df.index < s_date]
test_period = df.loc[(df.index >= s_date) & (df.index <= e_date)]


start_forecast = datetime.strptime(s_date, '%Y-%m-%d').date()
end_forecast = datetime.strptime(e_date, '%Y-%m-%d').date()


easter_start_2023 = '2023-04-01'
easter_end_2023 = '2023-04-22'

easter_weeks_test = df.loc[(df.index >= easter_start_2023) & (df.index <= easter_end_2023)]


from sklearn.metrics import mean_absolute_error, mean_squared_error
from math import sqrt
import numpy as np

def get_scores(actual, predicted):### Modelling Process
    # Calculate errors
    mae = mean_absolute_error(actual, predicted)
    mse = mean_squared_error(actual, predicted)
    rmse = sqrt(mse)
    percentage_diff = np.abs((actual - predicted) / actual) * 100

    # Calculate MAPE
    mape = percentage_diff.mean()

    # Calculate "Accuracy" Percentage
    accuracy_percentage = 100 - mape

    # Print metrics
    print(f"MAE: {mae}")
    print(f"MSE: {mse}")
    print(f"RMSE: {rmse}")
    print(f"MAPE: {mape}%")
    print(f"Direct 'Accuracy' Percentage: {accuracy_percentage}%")

    return pd.Series(data={'mae':mae, 'mse':mse, 'rmse':rmse, 'mape':mape, 'direct_accuracy':accuracy_percentage}, index=['mae', 'mse', 'rmse','mape','direct_accuracy'])


from statsmodels.tsa.arima.model import ARIMA

# ARIMA Model fitting
order = (1, 1, 1)  # Still (p, d, q)

model = ARIMA(train_period['UNIT_SALES'], order=order)
model_fit = model.fit()

# Forecasting
forecast = model_fit.predict(start=start_forecast,
                             end=end_forecast)

arima_forecast = forecast


arima_score = get_scores(test_period.squeeze(), arima_forecast)

MAE: 13878.577969407064
MSE: 333660932.83898854
RMSE: 18266.3880622029
MAPE: 13.029688220376284%
Direct 'Accuracy' Percentage: 86.97031177962371%


arima_easter_score = get_scores(easter_weeks_test.squeeze(), arima_forecast.loc[(arima_forecast.index >= easter_start_2023) & (arima_forecast.index <= easter_end_2023)])

MAE: 16694.848477498203
MSE: 332614641.0282627
RMSE: 18237.72576359955
MAPE: 15.726140646868652%
Direct 'Accuracy' Percentage: 84.27385935313134%


# Plotting the results
plt.figure(figsize=(10, 6))
plt.plot(easter_weeks_test['UNIT_SALES'], label='Actual Sales')
plt.plot(arima_forecast.loc[(arima_forecast.index >= easter_start_2023) & (arima_forecast.index <= easter_end_2023)], label='Forecast', color='red')
plt.xticks(easter_weeks_test.index)
plt.xlabel('Date')
plt.ylabel('UNIT_SALES')
plt.legend()
plt.show()


from statsmodels.tsa.statespace.sarimax import SARIMAX


# SARIMA Model fitting
order = (1, 1, 1)  # (p, d, q)
seasonal_order = (1, 1, 1, 52)  # (P, D, Q, s)

model = SARIMAX(train_period['UNIT_SALES'], order=order, seasonal_order=seasonal_order)
model_fit = model.fit()

# Forecasting
forecast = model_fit.predict(start=start_forecast.strftime('%Y-%m-%d'),
                             end=end_forecast.strftime('%Y-%m-%d'))
sarima_forecast = forecast


sarima_score = get_scores(test_period.squeeze(), sarima_forecast)

MAE: 10543.398089743641
MSE: 170655919.02185866
RMSE: 13063.533940777996
MAPE: 9.790936856842896%
Direct 'Accuracy' Percentage: 90.20906314315711%


sarima_easter_score = get_scores(easter_weeks_test.squeeze(), sarima_forecast.loc[(sarima_forecast.index >= easter_start_2023) & (sarima_forecast.index <= easter_end_2023)])

MAE: 19749.376548111453
MSE: 497937909.1387423
RMSE: 22314.522382043993
MAPE: 19.265721342075864%
Direct 'Accuracy' Percentage: 80.73427865792414%


# Plotting the results
plt.figure(figsize=(10, 6))
plt.plot(easter_weeks_test['UNIT_SALES'], label='Actual Sales')
plt.plot(sarima_forecast.loc[(sarima_forecast.index >= easter_start_2023) & (sarima_forecast.index <= easter_end_2023)], label='Forecast', color='red')
plt.xticks(easter_weeks_test.index)
plt.xlabel('Date')
plt.ylabel('UNIT_SALES')
plt.legend()
plt.show()


from prophet import Prophet

df_mod = train_period.reset_index()
df_mod['DATE'] = pd.to_datetime(df_mod['DATE'])
df_mod.columns = ['ds', 'y']

model = Prophet(changepoint_prior_scale=0.1,seasonality_prior_scale=0.5, weekly_seasonality=True)
model.fit(df_mod)

# Forecast the next 52 weeks (1 year)
future = model.make_future_dataframe(periods=53, freq='W-SAT')
forecast = model.predict(future)

prophet_forecast = forecast[['ds','yhat']]
prophet_forecast['ds'] = pd.to_datetime(prophet_forecast['ds'])
prophet_forecast.columns = ['DATE', 'yhat']
prophet_forecast.set_index('DATE', inplace=True)
prophet_forecast = prophet_forecast.asfreq('W-SAT')
prophet_forecast = prophet_forecast.loc[(prophet_forecast.index >= s_date) & (prophet_forecast.index <= e_date)]


prophet_score = get_scores(test_period.squeeze(), prophet_forecast.loc[(prophet_forecast.index >= s_date) & (prophet_forecast.index <= e_date)].squeeze())

MAE: 10354.142149221083
MSE: 178938113.5920231
RMSE: 13376.775156666987
MAPE: 9.585562452804588%
Direct 'Accuracy' Percentage: 90.4144375471954%


prophet_easter_score = get_scores(easter_weeks_test.squeeze(), prophet_forecast.loc[(prophet_forecast.index >= easter_start_2023) & (prophet_forecast.index <= easter_end_2023)].squeeze())

MAE: 13672.249500615759
MSE: 250724046.9852324
RMSE: 15834.268122816173
MAPE: 13.398864372818881%
Direct 'Accuracy' Percentage: 86.60113562718112%


plt.figure(figsize=(10, 6))
sns.lineplot(data=easter_weeks_test, x=easter_weeks_test.index, y='UNIT_SALES', label='Actual Sales')
sns.lineplot(data=prophet_forecast.loc[(prophet_forecast.index >= easter_start_2023) & (prophet_forecast.index <= easter_end_2023)], x=prophet_forecast.loc[(prophet_forecast.index >= easter_start_2023) & (prophet_forecast.index <= easter_end_2023)].index, y='yhat', label='Forecast', color='red')
plt.xticks(easter_weeks_test.index)
plt.xlabel('Date')
plt.ylabel('UNIT_SALES')
plt.legend()
plt.show()


from statsmodels.tsa.holtwinters import ExponentialSmoothing

# ExponentialSmoothing Model fitting
model = ExponentialSmoothing(train_period['UNIT_SALES'], trend='add', seasonal='add', seasonal_periods=52).fit(smoothing_level=0.6, smoothing_slope=0.2, smoothing_seasonal=0.2)

# Forecasting

forecast_periods = ((end_forecast - start_forecast).days // 7) + 1

exponential_forecast = model.forecast(forecast_periods)

forecast_dates = pd.date_range(start=start_forecast, periods=forecast_periods, freq='W-SAT')


exponential_smoothing_score = get_scores(test_period.squeeze(), exponential_forecast)

MAE: 11515.003782552438
MSE: 176593000.4153336
RMSE: 13288.829911445688
MAPE: 11.473932683680282%
Direct 'Accuracy' Percentage: 88.52606731631971%


exponential_smoothing_easter_score = get_scores(easter_weeks_test.squeeze(), exponential_forecast.loc[(exponential_forecast.index >= easter_start_2023) & (exponential_forecast.index <= easter_end_2023)])

MAE: 15194.8459825914
MSE: 339963988.23290277
RMSE: 18438.112382586856
MAPE: 15.957019954688251%
Direct 'Accuracy' Percentage: 84.04298004531175%


# Plotting the results

easter_start_2023 = pd.to_datetime(easter_start_2023)
easter_end_2023 = pd.to_datetime(easter_end_2023)

# Filter the dates
filtered_dates = forecast_dates[(forecast_dates >= easter_start_2023) & (forecast_dates <= easter_end_2023)]


plt.figure(figsize=(10, 6))
plt.plot(easter_weeks_test.index, easter_weeks_test['UNIT_SALES'], label='Actual Sales')
plt.plot(filtered_dates, exponential_forecast.loc[(exponential_forecast.index >= easter_start_2023) & (exponential_forecast.index <= easter_end_2023)], label='Forecast', color='red')
plt.xticks(easter_weeks_test.index)
plt.xlabel('Date')
plt.ylabel('UNIT_SALES')
plt.legend()
plt.show()


pd.options.display.float_format = '{:.2f}'.format
q2_scores = pd.DataFrame({'Arima':arima_score, 'Sarima':sarima_score, 'Prophet':prophet_score, 'Exponential Smoothing':exponential_smoothing_score}).T
print(q2_scores)

                           mae          mse     rmse  mape  direct_accuracy
Arima                 13878.58 333660932.84 18266.39 13.03            86.97
Sarima                10543.40 170655919.02 13063.53  9.79            90.21
Prophet               10354.14 178938113.59 13376.78  9.59            90.41
Exponential Smoothing 11515.00 176593000.42 13288.83 11.47            88.53


fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(18, 10))
fig.suptitle('Model Performance Metrics')

for i, column in enumerate(q2_scores.columns):
    row, col = divmod(i, 3)
    sns.barplot(ax=axes[row, col], x=q2_scores.index, y=q2_scores[column])
    axes[row, col].set_title(column)
    axes[row, col].set_ylabel('')

# Adjust layout
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()


q2_easter_scores = pd.DataFrame({'Arima':arima_easter_score, 'Sarima':sarima_easter_score, 'Prophet':prophet_easter_score, 'Exponential Smoothing':exponential_smoothing_easter_score}).T
print(q2_easter_scores)

                           mae          mse     rmse  mape  direct_accuracy
Arima                 16694.85 332614641.03 18237.73 15.73            84.27
Sarima                19749.38 497937909.14 22314.52 19.27            80.73
Prophet               13672.25 250724046.99 15834.27 13.40            86.60
Exponential Smoothing 15194.85 339963988.23 18438.11 15.96            84.04


fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(18, 10))
fig.suptitle('Model Performance Metrics')

for i, column in enumerate(q2_easter_scores.columns):
    row, col = divmod(i, 3)
    sns.barplot(ax=axes[row, col], x=q2_scores.index, y=q2_scores[column])
    axes[row, col].set_title(column)
    axes[row, col].set_ylabel('')

# Adjust layout
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()


df = avocado_regular_ssd_swire_small[['DATE','UNIT_SALES']].groupby(by="DATE", as_index=False).sum()
df['DATE'] = pd.to_datetime(df['DATE'])
df.set_index('DATE', inplace=True)
df = df.loc[df.index <= "2023-08-26"]
df = df.asfreq('W-SAT')


start_forecast = datetime.strptime("2024-03-23", '%Y-%m-%d').date()
end_forecast = datetime.strptime("2024-04-13", '%Y-%m-%d').date()


from prophet import Prophet

df_mod = df.reset_index()
df_mod['DATE'] = pd.to_datetime(df_mod['DATE'])
df_mod.columns = ['ds', 'y']

model = Prophet(changepoint_prior_scale=0.1,seasonality_prior_scale=0.5, weekly_seasonality=True)
model.fit(df_mod)

# Forecast the next 52 weeks (1 year)
future = model.make_future_dataframe(periods=53, freq='W-SAT')
forecast = model.predict(future)

prophet_forecast = forecast[['ds','yhat','yhat_upper','yhat_lower']]
prophet_forecast['ds'] = pd.to_datetime(prophet_forecast['ds'])
prophet_forecast.columns = ['DATE', 'yhat','yhat_upper','yhat_lower']
prophet_forecast.set_index('DATE', inplace=True)
prophet_forecast = prophet_forecast.asfreq('W-SAT')
prophet_forecast = prophet_forecast.loc[(prophet_forecast.index >= "2024-03-23") & (prophet_forecast.index <= "2024-04-13")]


prophet_forecast


plt.figure(figsize=(10, 6))
plt.plot(prophet_forecast.index, prophet_forecast['yhat'], label='Predicted', color='blue', marker='o')
plt.fill_between(prophet_forecast.index, prophet_forecast['yhat_lower'], prophet_forecast['yhat_upper'], color='gray', alpha=0.2, label='Confidence Interval')
plt.xticks(prophet_forecast.index)
plt.title('Forecast with Confidence Interval')
plt.xlabel('Date')
plt.ylabel('Forecasted UNIT_SALES')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()


# filtering for kiwano
kiwano = market_demand[(market_demand['ITEM'].str.contains('KIWANO', case=False, regex=True))]

# Filtering Swire + kiwano data
kiwano_swire = kiwano[kiwano['MANUFACTURER'] == 'SWIRE-CC']

# Filtering Swire +kiwano + Diet
kiwano_swire_diet = kiwano_swire[kiwano_swire['CALORIC_SEGMENT']  == 'DIET/LIGHT']


# Assuming 'DATE' column is in datetime format
kiwano_swire_diet['DATE'] = pd.to_datetime(kiwano_swire_diet['DATE'])
kiwano_swire_diet_sorted= kiwano_swire_diet.sort_values('DATE')

# Aggregate by Date and compute unit sales
kiwano_swire_diet_sorted_daily = kiwano_swire_diet.groupby('DATE')['UNIT_SALES'].sum().reset_index()


earliest_date = kiwano_swire_diet_sorted_daily['DATE'].min()

# Step 2: Calculate week numbers
kiwano_swire_diet_sorted_daily['WEEK_NUMBER'] = ((kiwano_swire_diet_sorted_daily['DATE'] - earliest_date).dt.days // 7) % 52

# Step 3: Aggregate the data for every week
weekly_data_kiwano = kiwano_swire_diet_sorted_daily.groupby('WEEK_NUMBER')['UNIT_SALES'].sum().reset_index()


import pandas as pd

# Assuming weekly_data is your DataFrame containing 'WEEK_NUMBER' and 'UNIT_SALES' columns

# Sort the DataFrame by 'WEEK_NUMBER' if it's not already sorted
weekly_data_kiwano = weekly_data_kiwano.sort_values(by='WEEK_NUMBER')

# Ensure that there are enough rows to compute the rolling sum
if len(weekly_data_kiwano) < 13:
    print("Insufficient data to compute rolling sum.")
else:
    # Initialize an empty list to store the results
    rolling_sales_sum = []

    # Iterate over the DataFrame to calculate the rolling sum for each 13-week interval
    for i in range(len(weekly_data_kiwano) - 12):
        # Calculate the sum of UNIT_SALES for the current 13-week interval
        sales_sum = weekly_data_kiwano.loc[i:i+12, 'UNIT_SALES'].sum()
        # Construct the 13-week interval string
        interval = f"{weekly_data_kiwano.iloc[i]['WEEK_NUMBER']}-{weekly_data_kiwano.iloc[i+12]['WEEK_NUMBER']}"
        # Append the result to the list
        rolling_sales_sum.append((sales_sum, interval))

    # Create a new DataFrame from the results
    rolling_sales_df = pd.DataFrame(rolling_sales_sum, columns=['SALES_13W_SUM', '13_WEEK_INTERVAL'])


import matplotlib.pyplot as plt

# Assuming rolling_sales_df contains the DataFrame with 'SALES_13W_SUM' and '13_WEEK_INTERVAL' columns

# Extract data for plotting
sales_sum = rolling_sales_df['SALES_13W_SUM']
interval = rolling_sales_df['13_WEEK_INTERVAL']

# Plot the rolling sum
plt.figure(figsize=(18, 6))
plt.plot(interval, sales_sum, marker='o', linestyle='-')
plt.xlabel('13-Week Interval')
plt.ylabel('Sales Sum')
plt.title('Rolling Sum of Sales over 13-Week Intervals')
plt.xticks(rotation=45)  # Rotate x-axis labels for better readability
plt.grid(True)

# Find the index of the maximum value
max_index = sales_sum.idxmax()
max_value = sales_sum[max_index]
max_interval = interval[max_index]

# Annotate the peak with a red marker and the corresponding 13-week interval
plt.annotate(f'Max: {max_value} ({max_interval})',
             xy=(max_interval, max_value),
             xytext=(70, -20),
             textcoords='offset points',
             arrowprops=dict(facecolor='red', arrowstyle='->'),
             fontsize=10,
             color='red',
             ha='center')

plt.tight_layout()
plt.show()


from prophet import Prophet
from prophet.plot import plot_plotly, plot_components_plotly
import pandas as pd



# Calculate the index for partitioning the data
partition_index = int(len(kiwano_swire_diet_sorted_daily) * 0.8)

# Split the data into train and test sets
train_data = kiwano_swire_diet_sorted_daily.iloc[:partition_index]
test_data = kiwano_swire_diet_sorted_daily.iloc[partition_index:]

# Rename columns as per Prophet's requirements
train_data = train_data.rename(columns={'DATE': 'ds', 'UNIT_SALES': 'y'})
test_data = test_data.rename(columns={'DATE': 'ds', 'UNIT_SALES': 'y'})

# Print the lengths of train and test dataframes
print("Length of Train Data:", len(train_data))
print("Length of Test Data:", len(test_data))

# Print the start and end dates of train and test dataframes
print("Start Date of Train Data:", train_data['ds'].min())
print("End Date of Train Data:", train_data['ds'].max())
print("Start Date of Test Data:", test_data['ds'].min())
print("End Date of Test Data:", test_data['ds'].max())

Length of Train Data: 118
Length of Test Data: 30
Start Date of Train Data: 2020-12-05 00:00:00
End Date of Train Data: 2023-03-04 00:00:00
Start Date of Test Data: 2023-03-11 00:00:00
End Date of Test Data: 2023-10-28 00:00:00


# Python
m = Prophet()
m.fit(train_data)
pred_df = m.predict(test_data)
pred_df.head()

12:27:59 - cmdstanpy - INFO - Chain [1] start processing
12:28:00 - cmdstanpy - INFO - Chain [1] done processing


fig1 = m.plot(pred_df)


from sklearn.metrics import mean_absolute_error, mean_squared_error
from math import sqrt
import numpy as np


def dividefn(num1, num2):
  try:
    return num1/num2
  except ZeroDivisionError:
    return 0


def mean_absolute_percentage_error(y_true, y_pred):
  y_true, y_pred = np.array(y_true), np.array(y_pred)
  out = [0] * len(y_true)
  for i in np.arange(len(y_true)):
    a = (y_pred[i] - y_true[i])*1
    b = y_true[i]
    out[i] = dividefn(a,b)
  return np.mean(np.abs(out))*100


def val_metrics_sales(inp_df):
  metric= []
  mape = np.round(mean_absolute_percentage_error(test_data['y'],pred_df['yhat']),3)
  mae = np.round(mean_absolute_error(test_data['y'], pred_df['yhat']),3)
  mse_error = np.round(mean_squared_error(test_data['y'],pred_df['yhat']),3)
  rmse = np.round(np.sqrt(mean_squared_error(test_data['y'],pred_df['yhat'])),3)
  accuracy = np.round((100 - abs(mape)),3)
  metric.append([mape,mae,mse_error,rmse,accuracy])
  metric = np.array(metric)
  return metric


Validation_metric = val_metrics_sales(pred_df)


print("MAPE:" , Validation_metric[0][0])
print("MAE:" , Validation_metric[0][1])
print("MSE_Error:" , Validation_metric[0][2])
print("RMSE:" , Validation_metric[0][3])
print("Accuracy:" , Validation_metric[0][4])

MAPE: 20.181
MAE: 3731.182
MSE_Error: 19613458.552
RMSE: 4428.708
Accuracy: 79.819


kiwano_swire_diet_sorted_daily.tail(1)


kiwano_future = kiwano_swire_diet_sorted_daily.copy()

#  Renaming columns for prophet requirements
kiwano_future.rename(columns = {'DATE':'ds','UNIT_SALES':'y'}, inplace = True)


m_2 = Prophet()
m_2.fit(kiwano_future)

# Create future dates for forecasting
future_kiwano = m_2.make_future_dataframe(periods=52,freq = 'W')

# Forecasting for 1 year
future_kiwano_forecast = m_2.predict(future_kiwano)

future_kiwano_forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail()

12:30:42 - cmdstanpy - INFO - Chain [1] start processing
12:30:42 - cmdstanpy - INFO - Chain [1] done processing


plot_plotly(m_2, future_kiwano_forecast)


# Filtering out prediction to get forecasted sales
forecast_future = future_kiwano_forecast[future_kiwano_forecast['ds'] > '2023-10-28']
forecast_future
# # No Negative Sales
forecast_future['yhat'] = forecast_future['yhat'].clip(lower=0)


forecast_future_date_sales = forecast_future[['ds','yhat']]


forecast_future_date_sales.head(5)


earliest_date = forecast_future_date_sales['ds'].min()

# Step 2: Calculate week numbers
forecast_future_date_sales['WEEK_NUMBER'] = ((forecast_future_date_sales['ds'] - earliest_date).dt.days // 7) % 52

# Step 3: Aggregate the data for every week
weekly_data_kiwano_forecasted = forecast_future_date_sales.groupby('WEEK_NUMBER')['yhat'].sum().reset_index()


import pandas as pd

# Assuming weekly_data is your DataFrame containing 'WEEK_NUMBER' and 'UNIT_SALES' columns

# Sort the DataFrame by 'WEEK_NUMBER' if it's not already sorted
weekly_data_kiwano_forecasted = weekly_data_kiwano_forecasted.sort_values(by='WEEK_NUMBER')

# Ensure that there are enough rows to compute the rolling sum
if len(weekly_data_kiwano_forecasted) < 13:
    print("Insufficient data to compute rolling sum.")
else:
    # Initialize an empty list to store the results
    rolling_sales_sum = []

    # Iterate over the DataFrame to calculate the rolling sum for each 13-week interval
    for i in range(len(weekly_data_kiwano_forecasted) - 12):
        # Calculate the sum of UNIT_SALES for the current 13-week interval
        sales_sum = weekly_data_kiwano_forecasted.loc[i:i+12, 'yhat'].sum()
        # Construct the 13-week interval string
        interval = f"{weekly_data_kiwano_forecasted.iloc[i]['WEEK_NUMBER']}-{weekly_data_kiwano_forecasted.iloc[i+12]['WEEK_NUMBER']}"
        # Append the result to the list
        rolling_sales_sum.append((sales_sum, interval))

    # Create a new DataFrame from the results
    rolling_sales_df = pd.DataFrame(rolling_sales_sum, columns=['SALES_13W_SUM', '13_WEEK_INTERVAL'])


import matplotlib.pyplot as plt

# Assuming rolling_sales_df contains the DataFrame with 'SALES_13W_SUM' and '13_WEEK_INTERVAL' columns

# Extract data for plotting
sales_sum = rolling_sales_df['SALES_13W_SUM']
interval = rolling_sales_df['13_WEEK_INTERVAL']

# Plot the rolling sum
plt.figure(figsize=(18, 6))
plt.plot(interval, sales_sum, marker='o', linestyle='-')
plt.xlabel('13-Week Interval')
plt.ylabel('Sales Sum')
plt.title('Rolling Sum of Sales over 13-Week Intervals')
plt.xticks(rotation=45)  # Rotate x-axis labels for better readability
plt.grid(True)

# Find the index of the maximum value
max_index = sales_sum.idxmax()
max_value = sales_sum[max_index]
max_interval = interval[max_index]

# Annotate the peak with a red marker and the corresponding 13-week interval
plt.annotate(f'Max: {max_value} ({max_interval})',
             xy=(max_interval, max_value),
             xytext=(70, -20),
             textcoords='offset points',
             arrowprops=dict(facecolor='red', arrowstyle='->'),
             fontsize=10,
             color='red',
             ha='center')

plt.tight_layout()
plt.show()


# Function to calculate the scores

from sklearn.metrics import mean_absolute_error, mean_squared_error
from math import sqrt
import numpy as np
def get_scores(actual, predicted):
# Calculate errors
    mae = mean_absolute_error(actual, predicted)
    mse = mean_squared_error(actual, predicted)
    rmse = sqrt(mse)
    #percentage_diff = (np.abs(actual - predicted).sum()) /( np.abs(actual).mean())
    
    # Combine actual and predicted values into one DataFrame
    mape_df = pd.DataFrame({'actual': actual.values, 'predicted': predicted.values})

# Calculate absolute percentage error
    mape_df['absolute_percentage_error'] = np.abs((mape_df['actual'] - mape_df['predicted']) / mape_df['actual']) * 100

# Calculate MAPE
    mape = mape_df['absolute_percentage_error'].mean()
    
# Calculate MAPE
    #mape = percentage_diff.mean()
# Calculate "Accuracy" Percentage
    accuracy_percentage = 100 - mape
# Print metrics
#     print(f"MAE: {mae}")
#     print(f"MSE: {mse}")
#     print(f"RMSE: {rmse}")
#     print(f"MAPE: {mape}%")
#     print(f"Direct 'Accuracy' Percentage: {accuracy_percentage}%")
    return pd.Series(data={'mae':mae, 'mse':mse, 'rmse':rmse, 'mape':mape, 'direct_accuracy':accuracy_percentage}, index=['mae', 'mse', 'rmse','mape','direct_accuracy'])


# Reading the zip_to_market data
zip_to_market = pd.read_csv(r"zip_to_market_unit_mapping.csv")
consumer_demographics = pd.read_csv(r"pivoted consumer data.csv")


# Dictionary mapping states to regions
state_regions = {
    'WA': 'North West',
    'OR': 'North West',
    'ID': 'North',
    'WY': 'North',
    'AZ': 'South',
    'NM': 'South',
    'NV': 'West',
    'CA': 'West',
    'UT': 'Central',
    'CO': 'Central',
    'NE': 'East',
    'SD': 'East',
    'KS': 'East'
}

# Function to get region for a state
def get_region(state):
    return state_regions.get(state, 'Unknown')

# Add a new column 'Region' based on 'State'
consumer_demographics['Region'] = consumer_demographics['State'].map(get_region)


# Selecting only Zip and Region column
consumer_demographics_zip_region = consumer_demographics[['Zip', 'Region']]

# Filtering the dataset to consider only North regions
consumer_demographics_north = consumer_demographics_zip_region[consumer_demographics_zip_region['Region'] == 'North']


# Perform inner join on 'Zip' and 'ZIP_CODE' columns
merged_data = pd.merge(consumer_demographics_north, zip_to_market, left_on='Zip', right_on='ZIP_CODE', how='inner')

# Remove duplicates in MARKET_KEY and Region
merged_data = merged_data.drop_duplicates(subset=['MARKET_KEY', 'Region'])

## Dropping Zip Column
merged_data.drop(columns='Zip', inplace=True)

# Dropping the Zip Code Column
merged_data.drop(columns='ZIP_CODE', inplace=True)


# Combining market_demand with zip_customer dataset
final_merged_data = pd.merge(market_demand, merged_data, on='MARKET_KEY', how='inner')

# Checking if there are null values
final_merged_data.isnull().sum().sum()

0


# Filtering the mulberry flavor from item
mulberry = final_merged_data[(final_merged_data['ITEM'].str.contains('MULBERRIES', case=False, regex=True))]
mulberry_swire = mulberry[mulberry['MANUFACTURER'] == 'SWIRE-CC']
mulberry_swire_diet = mulberry_swire[mulberry_swire['CALORIC_SEGMENT']  == 'DIET/LIGHT']
# Printing the shape of the new dataset which contains only mulberry
mulberry.shape

(12784, 11)


# Printing the shape
mulberry_swire_diet.shape

(1708, 11)


# Convert 'DATE' column to datetime format
mulberry_swire_diet['DATE'] = pd.to_datetime(mulberry_swire_diet['DATE'])

# Grouping by date column and computing unit sales to eliminate duplicates
weekly_data = mulberry_swire_diet.groupby('DATE')['UNIT_SALES'].sum().reset_index()

# Group by year and sum the UNIT_SALES
yearly_sales = weekly_data.groupby(weekly_data['DATE'].dt.year)['UNIT_SALES'].sum()

# Create a new DataFrame from the grouped data
yearly_sales_df = pd.DataFrame({'Year': yearly_sales.index, 'Unit_Sales': yearly_sales.values})

# Printing the dataframe
yearly_sales_df


import matplotlib.pyplot as plt

# Find the peak sales year and value
peak_sales_year = yearly_sales_df.loc[yearly_sales_df['Unit_Sales'].idxmax(), 'Year']
peak_sales_value = yearly_sales_df['Unit_Sales'].max()

# Plot
plt.figure(figsize=(10, 6))
plt.plot(yearly_sales_df['Year'], yearly_sales_df['Unit_Sales'], marker='o', linestyle='-')
plt.title('Yearly Unit Sales')
plt.xlabel('Year')
plt.ylabel('Unit Sales')
plt.grid(True)
plt.xticks(yearly_sales_df['Year'])  # Show all years on x-axis

# Mark the peak sales value
plt.annotate(f'Peak: {peak_sales_value} units', xy=(peak_sales_year, peak_sales_value),
             xytext=(peak_sales_year + 1, peak_sales_value + 1000),
             arrowprops=dict(facecolor='black', shrink=0.05))

plt.tight_layout()
plt.show()


# Converting  Date to datetime variable
mulberry_swire_diet['DATE'] = pd.to_datetime(mulberry_swire_diet['DATE'])

# Grouping by date to ensure only one record for each date
mulberry_swire_diet = mulberry_swire_diet.groupby('DATE')['UNIT_SALES'].sum().reset_index()


# Extracting mulberry flavor
mulberry = final_merged_data[(final_merged_data['ITEM'].str.contains('MULBERRIES', case=False, regex=True))]

# Filtering out the Swire sales
mulberry_swire = mulberry[mulberry['MANUFACTURER'] == 'SWIRE-CC']

# Filtering out swire+mulberry+diet
mulberry_swire_diet = mulberry_swire[mulberry_swire['CALORIC_SEGMENT']  == 'DIET/LIGHT']


# Copying the datset 
mulberry_swire_diet_forecast = mulberry_swire_diet.copy()


# Converting to date time
mulberry_swire_diet_forecast['DATE'] = pd.to_datetime(mulberry_swire_diet_forecast['DATE'])

# Grouping to ensure one record for each date
mulberry_swire_diet_forecast_daily = mulberry_swire_diet_forecast.groupby('DATE')['UNIT_SALES'].sum().reset_index()

# Sorting by date column
mulberry_swire_diet_forecast = mulberry_swire_diet_forecast_daily.sort_values(by='DATE')


# Renaming for Prophet model requirements
mulberry_swire_diet_forecast = mulberry_swire_diet_forecast.rename(columns={'DATE': 'ds', 'UNIT_SALES': 'y'})

# Extracting the last date
mulberry_swire_diet_forecast.tail(1)


# Creating an instance for the model
p_model1 = Prophet()

# Fitting the model
p_model1.fit(mulberry_swire_diet_forecast)

# Create future dates for forecasting for 1 year
future = p_model1.make_future_dataframe(periods=52,freq = 'W')

# Make predictions
forecast1 = p_model1.predict(future)

# Plot the forecast
# plot_plotly(p_model1, forecast1)


# Extracting the future sales for 1 year
forecast_future = forecast1[forecast1['ds'] > '2023-07-01']

# Printing to ensure if the filter is working right
forecast_future.head(2)


# Removing any negative sales
forecast_future['yhat'] = forecast_future['yhat'].clip(lower=0)


# Computing yearly sales
yearly_forecast_mulberry = forecast_future.groupby(forecast_future['ds'].dt.year)['yhat'].sum()

# Create a new DataFrame from the grouped data
yearly_forecast_mulberry_df = pd.DataFrame({'Year': yearly_forecast_mulberry.index, 'Unit_Sales': yearly_forecast_mulberry.values})

# Printing the data
yearly_forecast_mulberry_df


#### Extracing the mulberry flavor
mulberry = final_merged_data[(final_merged_data['ITEM'].str.contains('MULBERRIES', case=False, regex=True))]
mulberry_forecast = mulberry.copy()

# Converting to date time variable
mulberry_forecast['DATE'] = pd.to_datetime(mulberry_forecast['DATE'])

# No duplicate date rows
mulberry_forecast_daily = mulberry_forecast.groupby('DATE')['UNIT_SALES'].sum().reset_index()

# Sort according to the date
mulberry_forecast = mulberry_forecast_daily.sort_values(by='DATE')

# Rename for prophet requirements
mulberry_forecast = mulberry_forecast.rename(columns={'DATE': 'ds', 'UNIT_SALES': 'y'})

# Printing the last date 
mulberry_forecast.tail(5)


p_model2 = Prophet()
p_model2.fit(mulberry_forecast)

# Create future dates for forecasting
future = p_model2.make_future_dataframe(periods=52,freq = 'W')

# Make predictions
forecast2 = p_model2.predict(future)

# Plot the forecast
# plot_plotly(p_model2, forecast2)


# Extracting the forecated sales
forecast_future = forecast2[forecast2['ds'] > '2023-10-28']


forecast_future.tail(1)


forecast_future['yhat'] = forecast_future['yhat'].clip(lower=0)

yearly_forecast_mulberry = forecast_future.groupby(forecast_future['ds'].dt.year)['yhat'].sum()

# Create a new DataFrame from the grouped data
yearly_forecast_mulberry_df = pd.DataFrame({'Year': yearly_forecast_mulberry.index, 'Unit_Sales': yearly_forecast_mulberry.values})

yearly_forecast_mulberry_df


# Extracing mulberry flavor
mulberry = final_merged_data[(final_merged_data['ITEM'].str.contains('MULBERRIES', case=False, regex=True))]

# Filtering for diet 
mulberry_diet = mulberry[mulberry['CALORIC_SEGMENT']  == 'DIET/LIGHT']
mulberry_diet_forecast = mulberry_diet.copy()

# Converting to date time
mulberry_diet_forecast['DATE'] = pd.to_datetime(mulberry_diet_forecast['DATE'])

# Ensuring no duplicate records for date
mulberry_diet_forecast_daily = mulberry_diet_forecast.groupby('DATE')['UNIT_SALES'].sum().reset_index()

# Sorting the date
mulberry_diet_forecast_daily = mulberry_diet_forecast_daily.sort_values(by='DATE')

# Renaming for prophet requiremts
mulberry_diet_forecast = mulberry_diet_forecast_daily.rename(columns={'DATE': 'ds', 'UNIT_SALES': 'y'})

# Printing the last date
mulberry_diet_forecast.tail(3)


# Creating an instance for the model
p_model4 = Prophet()

# Fitting the model
p_model4.fit(mulberry_diet_forecast)

# Create future dates for forecasting
future = p_model4.make_future_dataframe(periods=52,freq = 'W')

# Make predictions
forecast3 = p_model4.predict(future)

# Plot the forecast
# plot_plotly(p_model4, forecast3)


# Extracting the future sales for 1 year
forecast_future = forecast3[forecast3['ds'] > '2023-10-28']


# Ensuring no negative unit sales
forecast_future['yhat'] = forecast_future['yhat'].clip(lower=0)

# Compute yearly sales
yearly_forecast_mulberry_diet = forecast_future.groupby(forecast_future['ds'].dt.year)['yhat'].sum()

# Create a new DataFrame from the grouped data
yearly_forecast_mulberry_diet_df = pd.DataFrame({'Year': yearly_forecast_mulberry_diet.index, 'Unit_Sales': yearly_forecast_mulberry_diet.values})
yearly_forecast_mulberry_diet_df


# Extracting mulberry flavor
mulberry = final_merged_data[(final_merged_data['ITEM'].str.contains('MULBERRIES', case=False, regex=True))]

# Filtering Swire sales
mulberry_swire = mulberry[mulberry['MANUFACTURER'] == 'SWIRE-CC']
mulberry_swire_forecast = mulberry_swire.copy()

# Converting to Date Time
mulberry_swire_forecast['DATE'] = pd.to_datetime(mulberry_swire_forecast['DATE'])

# Aggregating by date
mulberry_swire_forecast_daily = mulberry_swire_forecast.groupby('DATE')['UNIT_SALES'].sum().reset_index()

# Sorting by date
mulberry_swire_forecast = mulberry_swire_forecast_daily.sort_values(by='DATE')

# Renaming the columns for Prophet requirements
mulberry_swire_forecast = mulberry_swire_forecast.rename(columns={'DATE': 'ds', 'UNIT_SALES': 'y'})

# Printing the last date
mulberry_swire_forecast.tail(1)


# Creating an instance for the Prophet model
p_model3 = Prophet()

# Fitting the modle
p_model3.fit(mulberry_swire_forecast)

# Create future dates for forecasting
future = p_model3.make_future_dataframe(periods=52,freq = 'W')

# Make predictions
forecast2 = p_model3.predict(future)

# Plot the forecast
# plot_plotly(p_model3, forecast2)


# Future Sales
forecast_future = forecast2[forecast2['ds'] > '2023-07-01']

# No Negative Sales
forecast_future['yhat'] = forecast_future['yhat'].clip(lower=0)

# Year wise sales
yearly_forecast_mulberry_swire = forecast_future.groupby(forecast_future['ds'].dt.year)['yhat'].sum()

# Create a new DataFrame from the grouped data
yearly_forecast_mulberry_df = pd.DataFrame({'Year': yearly_forecast_mulberry_swire.index, 'Unit_Sales': yearly_forecast_mulberry_swire.values})
yearly_forecast_mulberry_df


mulberry = final_merged_data[(final_merged_data['ITEM'].str.contains('MULBERRIES', case=False, regex=True))]
mulberry_diet = mulberry[mulberry['CALORIC_SEGMENT']  == 'DIET/LIGHT']


mulberry_diet['CATEGORY'].value_counts()

CATEGORY
SPARKLING WATER       8399
ING ENHANCED WATER    3679
SSD                     18
Name: count, dtype: int64


mulberry_diet_sparkling = mulberry_diet[mulberry_diet['CATEGORY']  == 'SPARKLING WATER']


from sklearn.metrics import mean_absolute_error, mean_squared_error
from math import sqrt
import numpy as np
import pandas as pd
from prophet import Prophet
import matplotlib.pyplot as plt

# Assuming mulberry_diet_sparkling_validation is your DataFrame
mulberry_diet_sparkling_validation = mulberry_diet_sparkling.copy()

# Convert Date to datetime variable
mulberry_diet_sparkling_validation['DATE'] = pd.to_datetime(mulberry_diet_sparkling_validation['DATE'])

# Group by date to ensure only one record for each date
mulberry_diet_sparkling_validation = mulberry_diet_sparkling_validation.groupby('DATE')['UNIT_SALES'].sum().reset_index()

# Sort the data by date
mulberry_diet_sparkling_validation = mulberry_diet_sparkling_validation.sort_values(by='DATE')

# Calculate the index for partitioning the data
partition_index = int(len(mulberry_diet_sparkling_validation) * 0.8)

# Split the data into train and test sets
train_data = mulberry_diet_sparkling_validation.iloc[:partition_index]
test_data = mulberry_diet_sparkling_validation.iloc[partition_index:]

# Renaming for Prophet requirements
train_data = train_data.rename(columns={'DATE': 'ds', 'UNIT_SALES': 'y'})
test_data = test_data.rename(columns={'DATE': 'ds', 'UNIT_SALES': 'y'})

p_model = Prophet()
p_model.fit(train_data)

# Create future dates for forecasting
future = p_model.make_future_dataframe(periods=34, freq='W')

# Make predictions
forecast = p_model.predict(future)

# Extracting the date and unit sales from forecast
pred = forecast[['ds','yhat']]

# Plotting actual vs forecasted sales
plt.plot(train_data['ds'], train_data['y'], label='Actual Train Set', color='blue')
plt.plot(test_data['ds'], test_data['y'], label='Actual Test Set', color='green')
plt.plot(pred['ds'], pred['yhat'], label='Forecast', color='red')
plt.xlabel('Date')
plt.ylabel('Sales')
plt.title('Actual vs Forecast')
plt.legend()
plt.xticks(rotation=45)
plt.axhline(0, color='black', linewidth=0.5)
plt.show()

22:37:56 - cmdstanpy - INFO - Chain [1] start processing
22:37:56 - cmdstanpy - INFO - Chain [1] done processing


from sklearn.metrics import mean_absolute_error, mean_squared_error
from math import sqrt
import numpy as np
import pandas as pd
from prophet import Prophet
import matplotlib.pyplot as plt

# Assuming mulberry_diet_sparkling_validation is your DataFrame
mulberry_diet_sparkling_validation = mulberry_diet_sparkling.copy()

# Convert Date to datetime variable
mulberry_diet_sparkling_validation['DATE'] = pd.to_datetime(mulberry_diet_sparkling_validation['DATE'])

# Group by date to ensure only one record for each date
mulberry_diet_sparkling_validation = mulberry_diet_sparkling_validation.groupby('DATE')['UNIT_SALES'].sum().reset_index()

# Sort the data by date
mulberry_diet_sparkling_validation = mulberry_diet_sparkling_validation.sort_values(by='DATE')

# Calculate the index for partitioning the data
partition_index = int(len(mulberry_diet_sparkling_validation) * 0.8)

# Split the data into train and test sets
train_data = mulberry_diet_sparkling_validation.iloc[:partition_index]
test_data = mulberry_diet_sparkling_validation.iloc[partition_index:]

# Renaming for Prophet requirements
train_data = train_data.rename(columns={'DATE': 'ds', 'UNIT_SALES': 'y'})
test_data = test_data.rename(columns={'DATE': 'ds', 'UNIT_SALES': 'y'})

p_model = Prophet()
p_model.fit(train_data)

# Create future dates for forecasting
future = p_model.make_future_dataframe(periods=0, freq='W')

# Make predictions
forecast = p_model.predict(future)

# Extracting the date and unit sales from forecast
pred = forecast[['ds','yhat']]

get_scores(train_data['y'], pred['yhat'])

23:14:15 - cmdstanpy - INFO - Chain [1] start processing
23:14:15 - cmdstanpy - INFO - Chain [1] done processing

mae                  226.692952
mse                94751.882440
rmse                 307.817937
mape                   3.956043
direct_accuracy       96.043957
dtype: float64


forecast = p_model.predict(test_data)

pred = forecast[['ds','yhat']]
get_scores(test_data['y'], pred['yhat'])

mae                   654.519767
mse                618493.972874
rmse                  786.443878
mape                   13.653096
direct_accuracy        86.346904
dtype: float64


mulberry_diet_sparkling['DATE'] = pd.to_datetime(mulberry_diet_sparkling['DATE'])

# Ensuring no duplicate records for date
mulberry_diet_sparkling_daily = mulberry_diet_sparkling.groupby('DATE')['UNIT_SALES'].sum().reset_index()

# Sorting the date
mulberry_diet_sparkling_daily = mulberry_diet_sparkling_daily.sort_values(by='DATE')

# Renaming for prophet requiremts
mulberry_diet_sparkling_daily = mulberry_diet_sparkling_daily.rename(columns={'DATE': 'ds', 'UNIT_SALES': 'y'})

# Printing the last date
mulberry_diet_sparkling_daily.tail(1)


# Creating an instance for the model
p_model5 = Prophet()

# Fitting the model
p_model5.fit(mulberry_diet_sparkling_daily)

# Create future dates for forecasting
future5 = p_model5.make_future_dataframe(periods=52,freq = 'W')

# Make predictions
forecast5 = p_model5.predict(future5)

# Plot the forecast
plot_plotly(p_model5, forecast5)


# Future Sales
forecast5_future = forecast5[forecast5['ds'] > '2023-10-28']
forecast5_future
# # No Negative Sales
forecast5_future['yhat'] = forecast5_future['yhat'].clip(lower=0)

# # Year wise sales
yearly_forecast_mulberry_diet_sparkling = forecast5_future.groupby(forecast5_future['ds'].dt.year)['yhat'].sum()

# # Create a new DataFrame from the grouped data
yearly_forecast_mulberry_sparking_diet_df = pd.DataFrame({'Year': yearly_forecast_mulberry_diet_sparkling.index, 'Unit_Sales': yearly_forecast_mulberry_diet_sparkling.values})
yearly_forecast_mulberry_sparking_diet_df


mulberry = final_merged_data[(final_merged_data['ITEM'].str.contains('MULBERRIES', case=False, regex=True))]
mulberry_diet = mulberry[mulberry['CALORIC_SEGMENT']  == 'DIET/LIGHT']
mulberry_diet_sparkling = mulberry_diet[mulberry_diet['CATEGORY']  == 'SPARKLING WATER']


mulberry_diet_sparkling['MANUFACTURER'].value_counts(normalize=True)

MANUFACTURER
JOLLYS      0.571735
COCOS       0.151328
BEARS       0.141207
SWIRE-CC    0.135730
Name: proportion, dtype: float64


multi_jug_2l = market_demand[market_demand['PACKAGE']=='2L MULTI JUG']
diet_multi_jug_2l = multi_jug_2l[multi_jug_2l['CALORIC_SEGMENT']=='DIET/LIGHT']
swire_diet_multi_jug_2l = diet_multi_jug_2l[diet_multi_jug_2l['MANUFACTURER']=='SWIRE-CC']
swire_diet_multi_jug_2l_moonlit = swire_diet_multi_jug_2l[swire_diet_multi_jug_2l['BRAND'].str.contains('DIET MOONLIT', case=False, regex=True)]


df_6 = swire_diet_multi_jug_2l_moonlit[['DATE','UNIT_SALES']].groupby(by="DATE", as_index=False).sum()
df_6['DATE'] = pd.to_datetime(df_6['DATE'])
df_6.set_index('DATE', inplace=True)
df_6 = df_6.asfreq('W-SAT')


plt.plot(df_6)
plt.xlabel('Date')
plt.ylabel('UNIT_SALES')

Text(0, 0.5, 'UNIT_SALES')


# Train and Test periods
from datetime import datetime

s_date = '2023-02-25'
e_date = '2023-08-26'
df_6 = df_6.loc[df_6.index <= e_date]
train_period = df_6.loc[df_6.index < s_date]
test_period = df_6.loc[(df_6.index >= s_date) & (df_6.index <= e_date)]


start_forecast = datetime.strptime(s_date, '%Y-%m-%d').date()
end_forecast = datetime.strptime(e_date, '%Y-%m-%d').date()


from sklearn.metrics import mean_absolute_error, mean_squared_error
from math import sqrt
import numpy as np

def get_scores(actual, predicted):
    # Calculate errors
    mae = mean_absolute_error(actual, predicted)
    mse = mean_squared_error(actual, predicted)
    rmse = sqrt(mse)
    percentage_diff = np.abs((actual - predicted) / actual) * 100

    # Calculate MAPE
    mape = percentage_diff.mean()

    # Calculate "Accuracy" Percentage
    accuracy_percentage = 100 - mape

    # Print metrics
    print(f"MAE: {mae}")
    print(f"MSE: {mse}")
    print(f"RMSE: {rmse}")
    print(f"MAPE: {mape}%")
    print(f"Direct 'Accuracy' Percentage: {accuracy_percentage}%")

    return pd.Series(data={'mae':mae, 'mse':mse, 'rmse':rmse, 'mape':mape, 'direct_accuracy':accuracy_percentage}, index=['mae', 'mse', 'rmse','mape','direct_accuracy'])


from statsmodels.tsa.arima.model import ARIMA

# ARIMA Model fitting
order = (1, 1, 1)  # Still (p, d, q)

model = ARIMA(train_period['UNIT_SALES'], order=order)
model_fit = model.fit()

# Forecasting
forecast = model_fit.predict(start=start_forecast,
                             end=end_forecast)

arima_forecast = forecast


arima_score = get_scores(test_period.squeeze(), arima_forecast)

MAE: 1331.9101181487447
MSE: 2785766.362578469
RMSE: 1669.0615215079608
MAPE: 7.579412439517313%
Direct 'Accuracy' Percentage: 92.42058756048269%


# Plotting the results
plt.figure(figsize=(10, 6))
plt.plot(test_period['UNIT_SALES'], label='Actual Sales')
plt.plot(arima_forecast, label='Forecast', color='red')
plt.xlabel('Date')
plt.ylabel('UNIT_SALES')
plt.legend()
plt.show()


from statsmodels.tsa.statespace.sarimax import SARIMAX


# SARIMA Model fitting
order = (1, 1, 1)  # (p, d, q)
seasonal_order = (1, 1, 1, 52)  # (P, D, Q, s)

model = SARIMAX(train_period['UNIT_SALES'], order=order, seasonal_order=seasonal_order)
model_fit = model.fit()

# Forecasting
forecast = model_fit.predict(start=start_forecast.strftime('%Y-%m-%d'),
                             end=end_forecast.strftime('%Y-%m-%d'))


sarima_score = get_scores(test_period.squeeze(), sarima_forecast)

MAE: 2078.657764078491
MSE: 5833349.344575272
RMSE: 2415.2327723379526
MAPE: 12.472789935052829%
Direct 'Accuracy' Percentage: 87.52721006494717%


# Plotting the results
plt.figure(figsize=(10, 6))
plt.plot(test_period['UNIT_SALES'], label='Actual Sales')
plt.plot(sarima_forecast, label='Forecast', color='red')
plt.xlabel('Date')
plt.ylabel('UNIT_SALES')
plt.legend()
plt.show()


from prophet import Prophet

df_mod = train_period.reset_index()
df_mod['DATE'] = pd.to_datetime(df_mod['DATE'])
df_mod.columns = ['ds', 'y']

model = Prophet(changepoint_prior_scale=.1,seasonality_prior_scale=.1, weekly_seasonality=True)
model.fit(df_mod)

# Forecast the next 52 weeks (1 year)
future = model.make_future_dataframe(periods=53, freq='W-SAT')
forecast = model.predict(future)

prophet_forecast = forecast[['ds','yhat']]
prophet_forecast['ds'] = pd.to_datetime(prophet_forecast['ds'])
prophet_forecast.columns = ['DATE', 'yhat']
prophet_forecast.set_index('DATE', inplace=True)
prophet_forecast = prophet_forecast.asfreq('W-SAT')
prophet_forecast = prophet_forecast.loc[(prophet_forecast.index >= s_date) & (prophet_forecast.index <= e_date)]


prophet_score = get_scores(test_period.squeeze(), prophet_forecast.loc[(prophet_forecast.index >= s_date) & (prophet_forecast.index <= e_date)].squeeze())

MAE: 2573.16930571512
MSE: 9358981.759745052
RMSE: 3059.2452925100747
MAPE: 15.66648137188985%
Direct 'Accuracy' Percentage: 84.33351862811016%


# Plotting the results
plt.figure(figsize=(10, 6))
plt.plot(test_period['UNIT_SALES'], label='Actual Sales')
plt.plot(prophet_forecast, label='Forecast', color='red')
plt.xlabel('Date')
plt.ylabel('UNIT_SALES')
plt.legend()
plt.show()


from statsmodels.tsa.holtwinters import ExponentialSmoothing

# ExponentialSmoothing Model fitting
model = ExponentialSmoothing(train_period['UNIT_SALES'], trend='add', seasonal='add', seasonal_periods=52).fit(smoothing_level=0.01, smoothing_slope=0.01, smoothing_seasonal=0.1)

# Forecasting

forecast_periods = ((end_forecast - start_forecast).days // 7) + 1 

exponential_forecast = model.forecast(forecast_periods)

forecast_dates = pd.date_range(start=start_forecast, periods=forecast_periods, freq='W-SAT')


exponential_smoothing_score = get_scores(test_period.squeeze(), exponential_forecast)

MAE: 1672.5828724262556
MSE: 4512195.604211318
RMSE: 2124.192930082227
MAPE: 9.544250972143143%
Direct 'Accuracy' Percentage: 90.45574902785685%


# Plotting the results
plt.figure(figsize=(10, 6))
plt.plot(test_period['UNIT_SALES'], label='Actual Sales')
plt.plot(exponential_forecast, label='Forecast', color='red')
plt.xlabel('Date')
plt.ylabel('UNIT_SALES')
plt.legend()
plt.show()


pd.options.display.float_format = '{:.2f}'.format
q2_scores = pd.DataFrame({'Arima':arima_score, 'Sarima':sarima_score, 'Prophet':prophet_score, 'Exponential Smoothing':exponential_smoothing_score}).T
print(q2_scores)

                          mae        mse    rmse  mape  direct_accuracy
Arima                 1331.91 2785766.36 1669.06  7.58            92.42
Sarima                2078.66 5833349.34 2415.23 12.47            87.53
Prophet               2573.17 9358981.76 3059.25 15.67            84.33
Exponential Smoothing 1672.58 4512195.60 2124.19  9.54            90.46


fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(18, 10))
fig.suptitle('Model Performance Metrics')

for i, column in enumerate(q2_scores.columns):
    row, col = divmod(i, 3)
    sns.barplot(ax=axes[row, col], x=q2_scores.index, y=q2_scores[column])
    axes[row, col].set_title(column)
    axes[row, col].set_ylabel('')

# Adjust layout
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()


# grouping by date
df_6 = swire_diet_multi_jug_2l_moonlit[['DATE','UNIT_SALES']].groupby(by="DATE", as_index=False).sum()
df_6['DATE'] = pd.to_datetime(df_6['DATE'])
df_6.set_index('DATE', inplace=True)
df_6 = df_6.loc[df_6.index <= e_date]
df_6 = df_6.asfreq('W-SAT')


start_forecast = datetime.strptime("2024-09-02", '%Y-%m-%d').date()


from statsmodels.tsa.holtwinters import ExponentialSmoothing

# ExponentialSmoothing Model fitting
model = ExponentialSmoothing(df_6['UNIT_SALES'], trend='add', seasonal='add', seasonal_periods=52).fit(smoothing_level=0.01, smoothing_slope=0.01, smoothing_seasonal=0.1)

# Forecasting

forecast_periods = 26 

exponential_forecast = model.forecast(forecast_periods)

forecast_dates = pd.date_range(start=start_forecast, periods=forecast_periods, freq='W-SAT')


# dispalying the 6month forecast
exponential_forecast

2023-09-02   16170.76
2023-09-09   15186.27
2023-09-16   15677.02
2023-09-23   16999.64
2023-09-30   16727.67
2023-10-07   15485.82
2023-10-14   15160.90
2023-10-21   16290.53
2023-10-28   15278.34
2023-11-04   15843.64
2023-11-11   15831.58
2023-11-18   15807.30
2023-11-25   17151.41
2023-12-02   16644.21
2023-12-09   17270.94
2023-12-16   16978.95
2023-12-23   18534.64
2023-12-30   17596.79
2024-01-06   17048.72
2024-01-13   15438.83
2024-01-20   15804.36
2024-01-27   15687.73
2024-02-03   16612.23
2024-02-10   16418.08
2024-02-17   15646.95
2024-02-24   16181.15
Freq: W-SAT, dtype: float64


# Plotting the results
plt.figure(figsize=(10, 6))
plt.plot(exponential_forecast, label='Forecast', color='red')
plt.xlabel('Date')
plt.ylabel('UNIT_SALES')
plt.legend()
plt.show()


consumer_demographics_zip_region = consumer_demographics[['Zip', 'State','REGION']]

consumer_demographics_region = consumer_demographics_zip_region[(consumer_demographics_zip_region['REGION'] == 'Central') | (consumer_demographics_zip_region['REGION'] == 'North West') ]
consumer_demographics_region['REGION'].value_counts()

REGION
North West    990
Central       800
Name: count, dtype: int64


# (consumer_demographics_region['Zip'].value_counts()>1).sum()


zip_cons = pd.merge(consumer_demographics_region, zip_to_market, left_on='Zip', right_on='ZIP_CODE', how='inner')
zip_cons = zip_cons.drop_duplicates(subset=['MARKET_KEY', 'REGION'])


merged_data = pd.merge(market_demand, zip_cons, on='MARKET_KEY', how='inner')


merged_data = merged_data.drop(['Zip','ZIP_CODE'], axis=1)


merged_data['State'].value_counts(normalize=True)

State
CO    0.381271
WA    0.263675
OR    0.228770
UT    0.126284
Name: proportion, dtype: float64


# Group by 'State' and calculate summary statistics for 'UNIT_SALES'
unit_sales_summary_by_state = merged_data.groupby('State')['UNIT_SALES'].describe()

unit_sales_summary_by_state


print("unit sales:")
region_sales_stats = merged_data.groupby('REGION')['UNIT_SALES'].describe()
print(region_sales_stats)

# Group the data by 'REGION' and compute summary statistics for 'DOLLAR_SALES'
region_sales_summary = merged_data.groupby('REGION')['DOLLAR_SALES'].describe()

print("dollar sales:")
# Display the summary statistics
print(region_sales_summary)

unit sales:
                count        mean         std   min   25%   50%    75%      max
REGION                                                                         
Central     5870381.0  168.083194  438.258413  0.04  12.0  45.0  140.0  14773.0
North West  5695620.0  150.656967  343.386241  0.04  13.0  49.0  141.0  11100.0
dollar sales:
                count        mean          std   min    25%     50%     75%       max
REGION                                                                               
Central     5870381.0  575.904097  1464.522706  0.01  40.16  156.31  500.48  69472.26
North West  5695620.0  487.903103  1122.230124  0.01  44.68  159.52  458.65  46137.39


merged_data['DATE'] = pd.to_datetime(merged_data['DATE'])

earliest_date = merged_data['DATE'].min()

# Step 2: Calculate week numbers
merged_data['WEEK_NUMBER'] = ((merged_data['DATE'] - earliest_date).dt.days // 7) % 52

# Filter data for Central region
central = merged_data[merged_data['REGION'] == 'Central']

# Filter data for Northwest region
northwest = merged_data[merged_data['REGION'] == 'North West']


# Step 3: Aggregate the data for every week

northwest_fil_1 = northwest[(northwest['BRAND']=='GREETINGLE') & (northwest['CATEGORY']=='ING ENHANCED WATER')]
northwest_fil_1 = northwest_fil_1.groupby('DATE')['UNIT_SALES'].sum().reset_index()
# central_fil_1 = northwest[(northwest['BRAND']=='GREETINGLE') & (northwest['CATEGORY']=='ING ENHANCED WATER')]
# weekly_data = northwest_fil_1.groupby('WEEK_NUMBER')['UNIT_SALES'].sum().reset_index()
northwest_fil_1.tail(1)


from prophet import Prophet
from prophet.plot import plot_plotly, plot_components_plotly

# Initialize and fit the Prophet model
p_model = Prophet()

nw_model_data1 = northwest_fil_1.rename(columns={'DATE': 'ds', 'UNIT_SALES': 'y'})

p_model.fit(nw_model_data1)

# Create future dates for forecasting
future = p_model.make_future_dataframe(periods=52, freq='W')

# Make predictions
forecast = p_model.predict(future)

# Plot the forecast
plot_plotly(p_model, forecast)


# week wise sales for best 13 weeks
nw_fil1_forecast_1year = forecast[forecast['ds'] > '2023-10-28']
earliest_date = nw_fil1_forecast_1year['ds'].min()
nw_fil1_forecast_1year

# Step 2: Calculate week numbers
nw_fil1_forecast_1year['WEEK_NUMBER'] = ((nw_fil1_forecast_1year['ds'] - earliest_date).dt.days // 7) % 52
weekly_forecast_data = nw_fil1_forecast_1year.groupby('WEEK_NUMBER')['yhat'].sum().reset_index()


# Ensure that there are enough rows to compute the rolling sum
if len(weekly_forecast_data) < 13:
    print("Insufficient data to compute rolling sum.")
else:
    # Initialize an empty list to store the results
    rolling_sales_sum = []

    # Iterate over the DataFrame to calculate the rolling sum for each 13-week interval
    for i in range(len(weekly_forecast_data) - 12):
        # Calculate the sum of UNIT_SALES for the current 13-week interval
        sales_sum = weekly_forecast_data.loc[i:i+12, 'yhat'].sum()
        # Construct the 13-week interval string
        interval = f"{weekly_forecast_data.iloc[i]['WEEK_NUMBER']}-{weekly_forecast_data.iloc[i+12]['WEEK_NUMBER']}"
        # Append the result to the list
        rolling_sales_sum.append((sales_sum, interval))

    # Create a new DataFrame from the results
    rolling_sales_df = pd.DataFrame(rolling_sales_sum, columns=['SALES_13W_SUM', '13_WEEK_INTERVAL'])


plotting:

  Cell In[23], line 1
    plotting:
             ^
SyntaxError: invalid syntax


import matplotlib.pyplot as plt

# Assuming rolling_sales_df contains the DataFrame with 'SALES_13W_SUM' and '13_WEEK_INTERVAL' columns

# Extract data for plotting
sales_sum = rolling_sales_df['SALES_13W_SUM']
interval = rolling_sales_df['13_WEEK_INTERVAL']

# Plot the rolling sum
plt.figure(figsize=(18, 6))
plt.plot(interval, sales_sum, marker='o', linestyle='-')
plt.xlabel('13-Week Interval')
plt.ylabel('Sales Sum')
plt.title('Rolling Sum of Sales over 13-Week Intervals')
plt.xticks(rotation=45)  # Rotate x-axis labels for better readability
plt.grid(True)

# Find the index of the maximum value
max_index = sales_sum.idxmax()
max_value = sales_sum[max_index]
max_interval = interval[max_index]

# Annotate the peak with a red marker and the corresponding 13-week interval
plt.annotate(f'Max: {max_value} ({max_interval})',
             xy=(max_interval, max_value),
             xytext=(70, -20),
             textcoords='offset points',
             arrowprops=dict(facecolor='red', arrowstyle='->'),
             fontsize=10,
             color='red',
             ha='center')

plt.tight_layout()
plt.show()


total_fil1_nw = weekly_forecast_data['yhat'].sum()
total_13_week_sum = rolling_sales_df[rolling_sales_df['13_WEEK_INTERVAL']== '29.0-41.0'].sum()['SALES_13W_SUM']
percentage = total_13_week_sum*100/total_fil1_nw
print(round(percentage,2))

30.62


central_fil_1 = central[(central['BRAND']=='GREETINGLE') & (central['CATEGORY']=='ING ENHANCED WATER')]
central_fil_1 = central_fil_1.groupby('DATE')['UNIT_SALES'].sum().reset_index()
central_fil_1.tail(1)


from prophet import Prophet
from prophet.plot import plot_plotly, plot_components_plotly

# Initialize and fit the Prophet model
p_model_1 = Prophet()

cw_model_data1 = central_fil_1.rename(columns={'DATE': 'ds', 'UNIT_SALES': 'y'})

p_model_1.fit(cw_model_data1)

# Create future dates for forecasting
future_1 = p_model_1.make_future_dataframe(periods=52, freq='W')

# Make predictions
forecast_1 = p_model_1.predict(future_1)

# Plot the forecast
plot_plotly(p_model_1, forecast_1)


# week wise sales for best 13 weeks
cw_fil1_forecast_1year = forecast_1[forecast_1['ds'] > '2023-10-28']
earliest_date_1 = cw_fil1_forecast_1year['ds'].min()

# Step 2: Calculate week numbers
cw_fil1_forecast_1year['WEEK_NUMBER'] = ((cw_fil1_forecast_1year['ds'] - earliest_date_1).dt.days // 7) % 52
weekly_forecast_data_1 = cw_fil1_forecast_1year.groupby('WEEK_NUMBER')['yhat'].sum().reset_index()


# Ensure that there are enough rows to compute the rolling sum
if len(weekly_forecast_data_1) < 13:
    print("Insufficient data to compute rolling sum.")
else:
    # Initialize an empty list to store the results
    rolling_sales_sum = []

    # Iterate over the DataFrame to calculate the rolling sum for each 13-week interval
    for i in range(len(weekly_forecast_data_1) - 12):
        # Calculate the sum of UNIT_SALES for the current 13-week interval
        sales_sum = weekly_forecast_data_1.loc[i:i+12, 'yhat'].sum()
        # Construct the 13-week interval string
        interval = f"{weekly_forecast_data_1.iloc[i]['WEEK_NUMBER']}-{weekly_forecast_data_1.iloc[i+12]['WEEK_NUMBER']}"
        # Append the result to the list
        rolling_sales_sum.append((sales_sum, interval))

    # Create a new DataFrame from the results
    rolling_sales_df = pd.DataFrame(rolling_sales_sum, columns=['SALES_13W_SUM', '13_WEEK_INTERVAL'])


import matplotlib.pyplot as plt

# Assuming rolling_sales_df contains the DataFrame with 'SALES_13W_SUM' and '13_WEEK_INTERVAL' columns

# Extract data for plotting
sales_sum = rolling_sales_df['SALES_13W_SUM']
interval = rolling_sales_df['13_WEEK_INTERVAL']

# Plot the rolling sum
plt.figure(figsize=(18, 6))
plt.plot(interval, sales_sum, marker='o', linestyle='-')
plt.xlabel('13-Week Interval')
plt.ylabel('Sales Sum')
plt.title('Rolling Sum of Sales over 13-Week Intervals')
plt.xticks(rotation=45)  # Rotate x-axis labels for better readability
plt.grid(True)

# Find the index of the maximum value
max_index = sales_sum.idxmax()
max_value = sales_sum[max_index]
max_interval = interval[max_index]

# Annotate the peak with a red marker and the corresponding 13-week interval
plt.annotate(f'Max: {max_value} ({max_interval})',
             xy=(max_interval, max_value),
             xytext=(70, -20),
             textcoords='offset points',
             arrowprops=dict(facecolor='red', arrowstyle='->'),
             fontsize=10,
             color='red',
             ha='center')

plt.tight_layout()
plt.show()


total_fil1_central = weekly_forecast_data_1['yhat'].sum()
total_13_week_sum_central = rolling_sales_df[rolling_sales_df['13_WEEK_INTERVAL']== '25.0-37.0'].sum()['SALES_13W_SUM']
percentage = total_13_week_sum_central*100/total_fil1_central
print(round(percentage,2))

31.24


northwest = merged_data[merged_data['REGION'] == 'North West']
woodsy_yellow = northwest[(northwest['ITEM'].str.contains('WOODSY  YELLOW ', case=False, regex=True))]
woodsy_yellow_regular = woodsy_yellow[woodsy_yellow['CALORIC_SEGMENT'] == 'REGULAR']
woodsy_yellow_regular_enhanced_water = woodsy_yellow_regular[woodsy_yellow_regular['CATEGORY'] == 'ING ENHANCED WATER']


woodsy_yellow_regular_enhanced_water['MANUFACTURER'].value_counts()

MANUFACTURER
COCOS    12285
Name: count, dtype: int64


northwest_woodsy = woodsy_yellow_regular_enhanced_water.groupby('DATE')['UNIT_SALES'].sum().reset_index()
northwest_woodsy.tail(1)


from prophet import Prophet
from prophet.plot import plot_plotly, plot_components_plotly

# Initialize and fit the Prophet model
p_model_3 = Prophet()

northwest_woodsy = northwest_woodsy.rename(columns={'DATE': 'ds', 'UNIT_SALES': 'y'})

p_model_3.fit(northwest_woodsy)

# Create future dates for forecasting
future_3 = p_model_3.make_future_dataframe(periods=52, freq='W')

# Make predictions
forecast_3 = p_model_3.predict(future_3)

# Plot the forecast
plot_plotly(p_model_3, forecast_3)


# week wise sales for best 13 weeks
northwest_woodsy_forecast_1year = forecast_3[forecast_3['ds'] > '2023-10-28']
earliest_date_4 = northwest_woodsy_forecast_1year['ds'].min()

# Step 2: Calculate week numbers
northwest_woodsy_forecast_1year['WEEK_NUMBER'] = ((northwest_woodsy_forecast_1year['ds'] - earliest_date_4).dt.days // 7) % 52
weekly_forecast_data_northwest_woodsy = northwest_woodsy_forecast_1year.groupby('WEEK_NUMBER')['yhat'].sum().reset_index()


# Ensure that there are enough rows to compute the rolling sum
if len(weekly_forecast_data_northwest_woodsy) < 13:
    print("Insufficient data to compute rolling sum.")
else:
    # Initialize an empty list to store the results
    rolling_sales_sum = []

    # Iterate over the DataFrame to calculate the rolling sum for each 13-week interval
    for i in range(len(weekly_forecast_data_northwest_woodsy) - 12):
        # Calculate the sum of UNIT_SALES for the current 13-week interval
        sales_sum = weekly_forecast_data_northwest_woodsy.loc[i:i+12, 'yhat'].sum()
        # Construct the 13-week interval string
        interval = f"{weekly_forecast_data_northwest_woodsy.iloc[i]['WEEK_NUMBER']}-{weekly_forecast_data_northwest_woodsy.iloc[i+12]['WEEK_NUMBER']}"
        # Append the result to the list
        rolling_sales_sum.append((sales_sum, interval))

    # Create a new DataFrame from the results
    rolling_sales_df = pd.DataFrame(rolling_sales_sum, columns=['SALES_13W_SUM', '13_WEEK_INTERVAL'])


import matplotlib.pyplot as plt

# Assuming rolling_sales_df contains the DataFrame with 'SALES_13W_SUM' and '13_WEEK_INTERVAL' columns

# Extract data for plotting
sales_sum = rolling_sales_df['SALES_13W_SUM']
interval = rolling_sales_df['13_WEEK_INTERVAL']

# Plot the rolling sum
plt.figure(figsize=(18, 6))
plt.plot(interval, sales_sum, marker='o', linestyle='-')
plt.xlabel('13-Week Interval')
plt.ylabel('Sales Sum')
plt.title('Rolling Sum of Sales over 13-Week Intervals')
plt.xticks(rotation=45)  # Rotate x-axis labels for better readability
plt.grid(True)

# Find the index of the maximum value
max_index = sales_sum.idxmax()
max_value = sales_sum[max_index]
max_interval = interval[max_index]

# Annotate the peak with a red marker and the corresponding 13-week interval
plt.annotate(f'Max: {max_value} ({max_interval})',
             xy=(max_interval, max_value),
             xytext=(70, -20),
             textcoords='offset points',
             arrowprops=dict(facecolor='red', arrowstyle='->'),
             fontsize=10,
             color='red',
             ha='center')

plt.tight_layout()
plt.show()


total_fil1_northwest_woodsy = weekly_forecast_data_northwest_woodsy['yhat'].sum()
total_13_week_sum_northwest_woodsy = rolling_sales_df[rolling_sales_df['13_WEEK_INTERVAL']== '30.0-42.0'].sum()['SALES_13W_SUM']
percentage = total_13_week_sum_northwest_woodsy*100/total_fil1_northwest_woodsy
print(round(percentage,2))

30.98


central = merged_data[merged_data['REGION'] == 'Central']
woodsy_yellow_central = central[(central['ITEM'].str.contains('WOODSY  YELLOW ', case=False, regex=True))]
woodsy_yellow_central_regular = woodsy_yellow_central[woodsy_yellow_central['CALORIC_SEGMENT'] == 'REGULAR']
woodsy_yellow_central_regular_enhanced_water = woodsy_yellow_central_regular[woodsy_yellow_central_regular['CATEGORY'] == 'ING ENHANCED WATER']


woodsy_yellow_central_regular_enhanced_water['MANUFACTURER'].value_counts()

MANUFACTURER
COCOS    11438
Name: count, dtype: int64


central_woodsy = woodsy_yellow_central_regular_enhanced_water.groupby('DATE')['UNIT_SALES'].sum().reset_index()
central_woodsy.tail(1)


from prophet import Prophet
from prophet.plot import plot_plotly, plot_components_plotly

# Initialize and fit the Prophet model
p_model_4 = Prophet()

central_woodsy = central_woodsy.rename(columns={'DATE': 'ds', 'UNIT_SALES': 'y'})

p_model_4.fit(central_woodsy)

# Create future dates for forecasting
future_5 = p_model_4.make_future_dataframe(periods=52, freq='W')

# Make predictions
forecast_5 = p_model_4.predict(future_5)

# Plot the forecast
plot_plotly(p_model_4, forecast_5)


# week wise sales for best 13 weeks
central_woodsy_forecast_1year = forecast_5[forecast_5['ds'] > '2023-10-28']
earliest_date_6 = central_woodsy_forecast_1year['ds'].min()

# Step 2: Calculate week numbers
central_woodsy_forecast_1year['WEEK_NUMBER'] = ((central_woodsy_forecast_1year['ds'] - earliest_date_6).dt.days // 7) % 52
weekly_forecast_data_central_woodsy = central_woodsy_forecast_1year.groupby('WEEK_NUMBER')['yhat'].sum().reset_index()


# Ensure that there are enough rows to compute the rolling sum
if len(weekly_forecast_data_central_woodsy) < 13:
    print("Insufficient data to compute rolling sum.")
else:
    # Initialize an empty list to store the results
    rolling_sales_sum = []

    # Iterate over the DataFrame to calculate the rolling sum for each 13-week interval
    for i in range(len(weekly_forecast_data_central_woodsy) - 12):
        # Calculate the sum of UNIT_SALES for the current 13-week interval
        sales_sum = weekly_forecast_data_central_woodsy.loc[i:i+12, 'yhat'].sum()
        # Construct the 13-week interval string
        interval = f"{weekly_forecast_data_central_woodsy.iloc[i]['WEEK_NUMBER']}-{weekly_forecast_data_central_woodsy.iloc[i+12]['WEEK_NUMBER']}"
        # Append the result to the list
        rolling_sales_sum.append((sales_sum, interval))

    # Create a new DataFrame from the results
    rolling_sales_df = pd.DataFrame(rolling_sales_sum, columns=['SALES_13W_SUM', '13_WEEK_INTERVAL'])


import matplotlib.pyplot as plt

# Assuming rolling_sales_df contains the DataFrame with 'SALES_13W_SUM' and '13_WEEK_INTERVAL' columns

# Extract data for plotting
sales_sum = rolling_sales_df['SALES_13W_SUM']
interval = rolling_sales_df['13_WEEK_INTERVAL']

# Plot the rolling sum
plt.figure(figsize=(18, 6))
plt.plot(interval, sales_sum, marker='o', linestyle='-')
plt.xlabel('13-Week Interval')
plt.ylabel('Sales Sum')
plt.title('Rolling Sum of Sales over 13-Week Intervals')
plt.xticks(rotation=45)  # Rotate x-axis labels for better readability
plt.grid(True)

# Find the index of the maximum value
max_index = sales_sum.idxmax()
max_value = sales_sum[max_index]
max_interval = interval[max_index]

# Annotate the peak with a red marker and the corresponding 13-week interval
plt.annotate(f'Max: {max_value} ({max_interval})',
             xy=(max_interval, max_value),
             xytext=(70, -20),
             textcoords='offset points',
             arrowprops=dict(facecolor='red', arrowstyle='->'),
             fontsize=10,
             color='red',
             ha='center')

plt.tight_layout()
plt.show()


total_fil1_central_woodsy = weekly_forecast_data_central_woodsy['yhat'].sum()
total_13_week_sum_northwest_central = rolling_sales_df[rolling_sales_df['13_WEEK_INTERVAL']== '32.0-44.0'].sum()['SALES_13W_SUM']
percentage = total_13_week_sum_northwest_central*100/total_fil1_central_woodsy
print(round(percentage,2))

31.14


import pandas as pd

data = {
    'Region': ['NorthWest', 'Central', 'NorthWest', 'Central'],
    'Filter': ['Brand  + Category', 'Brand  + Category ', 'Flavor + Caloric + Category ', 'Flavor  + Caloric + Category '],
    'Best 13 Weeks (Forecasted)': ['29 - 41', '25 - 37', '30 - 42', '32 - 44'],
    'Absolute Value (Unit Sales)': [300731, 272066, 173096, 108107],
    'Absolute Value Percentage for total weeks': ['30.62%', '31.24%', '30.98%', '31.14%']
}

df = pd.DataFrame(data)

df


southern_keys = [
    259, 260, 524, 272, 531, 277, 534, 535, 793, 794, 799, 801, 294, 807, 44,
    557, 310, 825, 316, 831, 323, 585, 849, 862, 866, 867, 358, 365, 622, 367,
    882, 891, 893, 895, 900, 913, 409, 673, 674, 939, 179, 949, 951, 952, 953,
    197, 967, 205, 210, 979, 212, 980, 981, 216, 729, 220, 733, 478, 739, 743,
    746, 748, 238, 754, 759, 766
]

# Filter the DataFrame
southern = market_demand[market_demand['MARKET_KEY'].isin(southern_keys)]


peppy = southern[(southern['BRAND'].str.contains('peppy', case=False, regex=True))]
peppy_regular = peppy[(peppy['CALORIC_SEGMENT'].str.contains('Regular', case=False, regex=True))]
peppy_regular_ssd = peppy_regular[(peppy_regular['CALORIC_SEGMENT'].str.contains('Regular', case=False, regex=True))]
peppy_regular_ssd_multi_jug = peppy_regular_ssd[(peppy_regular_ssd['PACKAGE'].str.contains('Multi Jug', case=False, regex=True))]


df = peppy_regular_ssd_multi_jug[['DATE','UNIT_SALES']].groupby(by="DATE", as_index=False).sum()
df['DATE'] = pd.to_datetime(df['DATE'])
df.set_index('DATE', inplace=True)
df = df.asfreq('W-SAT')


plt.plot(df)

[<matplotlib.lines.Line2D at 0x290d3368b90>]


# Train and Test periods
from datetime import datetime

s_date = '2023-06-03'
e_date = '2023-08-26'
df = df.loc[df.index <= e_date]
train_period = df.loc[df.index < s_date]
test_period = df.loc[(df.index >= s_date) & (df.index <= e_date)]


start_forecast = datetime.strptime(s_date, '%Y-%m-%d').date()
end_forecast = datetime.strptime(e_date, '%Y-%m-%d').date()


from sklearn.metrics import mean_absolute_error, mean_squared_error
from math import sqrt
import numpy as np

def get_scores(actual, predicted):
    # Calculate errors
    mae = mean_absolute_error(actual, predicted)
    mse = mean_squared_error(actual, predicted)
    rmse = sqrt(mse)
    percentage_diff = np.abs((actual - predicted) / actual) * 100

    # Calculate MAPE
    mape = percentage_diff.mean()

    # Calculate "Accuracy" Percentage
    accuracy_percentage = 100 - mape

    # Print metrics
    print(f"MAE: {mae}")
    print(f"MSE: {mse}")
    print(f"RMSE: {rmse}")
    print(f"MAPE: {mape}%")
    print(f"Direct 'Accuracy' Percentage: {accuracy_percentage}%")

    return pd.Series(data={'mae':mae, 'mse':mse, 'rmse':rmse, 'mape':mape, 'direct_accuracy':accuracy_percentage}, index=['mae', 'mse', 'rmse','mape','direct_accuracy'])


from statsmodels.tsa.arima.model import ARIMA

# ARIMA Model fitting
order = (2, 1, 2)  # Still (p, d, q)

model = ARIMA(train_period['UNIT_SALES'], order=order)
model_fit = model.fit()

# Forecasting
forecast = model_fit.predict(start=start_forecast,
                             end=end_forecast)

arima_forecast = forecast


arima_score = get_scores(test_period.squeeze(), arima_forecast)

MAE: 10610.921900092897
MSE: 133607735.47205538
RMSE: 11558.881237907732
MAPE: 3.682049809516563%
Direct 'Accuracy' Percentage: 96.31795019048344%


# Plotting the results
plt.figure(figsize=(10, 6))
plt.plot(test_period['UNIT_SALES'], label='Actual Sales')
plt.plot(arima_forecast, label='Forecast', color='red')
plt.xlabel('Date')
plt.ylabel('UNIT_SALES')
plt.legend()
plt.show()


from statsmodels.tsa.statespace.sarimax import SARIMAX


# SARIMA Model fitting
order = (1, 1, 1)  # (p, d, q)
seasonal_order = (1, 1, 1, 52)  # (P, D, Q, s)

model = SARIMAX(train_period['UNIT_SALES'], order=order, seasonal_order=seasonal_order)
model_fit = model.fit()

# Forecasting
forecast = model_fit.predict(start=start_forecast.strftime('%Y-%m-%d'),
                             end=end_forecast.strftime('%Y-%m-%d'))


sarima_score = get_scores(test_period.squeeze(), sarima_forecast)

MAE: 24348.214490805636
MSE: 801042551.4448053
RMSE: 28302.695126874496
MAPE: 8.420642444814801%
Direct 'Accuracy' Percentage: 91.5793575551852%


# Plotting the results
plt.figure(figsize=(10, 6))
plt.plot(test_period['UNIT_SALES'], label='Actual Sales')
plt.plot(sarima_forecast, label='Forecast', color='red')
plt.xlabel('Date')
plt.ylabel('UNIT_SALES')
plt.legend()
plt.show()


from prophet import Prophet

df_mod = train_period.reset_index()
df_mod['DATE'] = pd.to_datetime(df_mod['DATE'])
df_mod.columns = ['ds', 'y']

model = Prophet(changepoint_prior_scale=0.1,seasonality_prior_scale=0.5, weekly_seasonality=True)
model.fit(df_mod)

# Forecast the next 52 weeks (1 year)
future = model.make_future_dataframe(periods=53, freq='W-SAT')
forecast = model.predict(future)

prophet_forecast = forecast[['ds','yhat']]
prophet_forecast['ds'] = pd.to_datetime(prophet_forecast['ds'])
prophet_forecast.columns = ['DATE', 'yhat']
prophet_forecast.set_index('DATE', inplace=True)
prophet_forecast = prophet_forecast.asfreq('W-SAT')
prophet_forecast = prophet_forecast.loc[(prophet_forecast.index >= s_date) & (prophet_forecast.index <= e_date)]


prophet_score = get_scores(test_period.squeeze(), prophet_forecast.loc[(prophet_forecast.index >= s_date) & (prophet_forecast.index <= e_date)].squeeze())

MAE: 14082.66495499823
MSE: 257837638.2131654
RMSE: 16057.32350714668
MAPE: 4.88229507297955%
Direct 'Accuracy' Percentage: 95.11770492702045%


# Plotting the results
plt.figure(figsize=(10, 6))
plt.plot(test_period['UNIT_SALES'], label='Actual Sales')
plt.plot(prophet_forecast, label='Forecast', color='red')
plt.xlabel('Date')
plt.ylabel('UNIT_SALES')
plt.legend()
plt.show()


from statsmodels.tsa.holtwinters import ExponentialSmoothing

# ExponentialSmoothing Model fitting
model = ExponentialSmoothing(train_period['UNIT_SALES'], trend='add', seasonal='add', seasonal_periods=52).fit(smoothing_level=0.6, smoothing_slope=0.2, smoothing_seasonal=0.2)

# Forecasting

forecast_periods = ((end_forecast - start_forecast).days // 7) + 1

exponential_forecast = model.forecast(forecast_periods)

forecast_dates = pd.date_range(start=start_forecast, periods=forecast_periods, freq='W-SAT')


exponential_smoothing_score = get_scores(test_period.squeeze(), exponential_forecast)

MAE: 23907.5251431621
MSE: 744489010.3472316
RMSE: 27285.325916089616
MAPE: 8.281600968319658%
Direct 'Accuracy' Percentage: 91.71839903168035%


# Plotting the results
plt.figure(figsize=(10, 6))
plt.plot(test_period['UNIT_SALES'], label='Actual Sales')
plt.plot(exponential_forecast, label='Forecast', color='red')
plt.xlabel('Date')
plt.ylabel('UNIT_SALES')
plt.legend()
plt.show()


pd.options.display.float_format = '{:.2f}'.format
q2_scores = pd.DataFrame({'Arima':arima_score, 'Sarima':sarima_score, 'Prophet':prophet_score, 'Exponential Smoothing':exponential_smoothing_score}).T
print(q2_scores)

                           mae          mse     rmse  mape  direct_accuracy
Arima                 10610.92 133607735.47 11558.88  3.68            96.32
Sarima                24348.21 801042551.44 28302.70  8.42            91.58
Prophet               14082.66 257837638.21 16057.32  4.88            95.12
Exponential Smoothing 23907.53 744489010.35 27285.33  8.28            91.72


fig, axes = plt.subplots(nrows=2, ncols=3, figsize=(18, 10))
fig.suptitle('Model Performance Metrics')

for i, column in enumerate(q2_scores.columns):
    row, col = divmod(i, 3)
    sns.barplot(ax=axes[row, col], x=q2_scores.index, y=q2_scores[column])
    axes[row, col].set_title(column)
    axes[row, col].set_ylabel('')  # Removing the y-label for cleanliness

# Adjust layout
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()


df = peppy_regular_ssd_multi_jug[['DATE','UNIT_SALES']].groupby(by="DATE", as_index=False).sum()
df['DATE'] = pd.to_datetime(df['DATE'])
df.set_index('DATE', inplace=True)
df = df.loc[df.index <= "2023-08-26"]
df = df.asfreq('W-SAT')


from prophet import Prophet

df_mod = df.reset_index()
df_mod['DATE'] = pd.to_datetime(df_mod['DATE'])
df_mod.columns = ['ds', 'y']

model = Prophet(changepoint_prior_scale=0.1,seasonality_prior_scale=0.5, weekly_seasonality=True)
model.fit(df_mod)

future = model.make_future_dataframe(periods=13, freq='W-SAT')
forecast = model.predict(future)

prophet_forecast = forecast[['ds','yhat','yhat_upper','yhat_lower']]
prophet_forecast['ds'] = pd.to_datetime(prophet_forecast['ds'])
prophet_forecast.columns = ['DATE', 'yhat','yhat_upper','yhat_lower']
prophet_forecast.set_index('DATE', inplace=True)
prophet_forecast = prophet_forecast.asfreq('W-SAT')
prophet_forecast = prophet_forecast.loc[(prophet_forecast.index >= "2023-09-02")]

20:33:04 - cmdstanpy - INFO - Chain [1] start processing
20:33:04 - cmdstanpy - INFO - Chain [1] done processing
C:\Users\Michael Mendoza\AppData\Local\Temp\ipykernel_19324\963624670.py:14: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  prophet_forecast['ds'] = pd.to_datetime(prophet_forecast['ds'])


prophet_forecast


plt.figure(figsize=(10, 6))
plt.plot(prophet_forecast.index, prophet_forecast['yhat'], label='Predicted', color='blue', marker='o')
plt.fill_between(prophet_forecast.index, prophet_forecast['yhat_lower'], prophet_forecast['yhat_upper'], color='gray', alpha=0.2, label='Confidence Interval')
plt.title('Forecast with Confidence Interval')
plt.xlabel('Date')
plt.ylabel('Forecast Value')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

	ds	trend	yhat_lower	yhat_upper	trend_lower	trend_upper	additive_terms	additive_terms_lower	additive_terms_upper	yearly	yearly_lower	yearly_upper	multiplicative_terms	multiplicative_terms_lower	multiplicative_terms_upper	yhat	WEEK_NUMBER
145	2023-12-03	752.538260	476.851987	886.220095	750.513944	754.733644	-61.95856	-61.95856	-61.95856	-61.95856	-61.95856	-61.95856	0.0	0.0	0.0	690.579701	0
146	2023-12-10	749.318607	435.234791	850.450446	746.771570	752.273543	-89.59627	-89.59627	-89.59627	-89.59627	-89.59627	-89.59627	0.0	0.0	0.0	659.722337	1

	SSD	ENERGY
count	1.291228e+07	5.932087e+06
mean	2.088707e+02	2.020072e+02
std	1.026061e+03	8.425682e+02
min	4.000000e-02	4.000000e-02
25%	1.500000e+01	7.000000e+00
50%	5.300000e+01	3.300000e+01
75%	1.560000e+02	1.370000e+02
max	9.677600e+04	7.165700e+04

	DIET SMASH	VENOMOUS BLAST
count	17483.000000	51756.000000
mean	27.975093	69.975983
std	30.263642	246.676148
min	1.000000	1.000000
25%	9.000000	6.000000
50%	21.000000	15.000000
75%	38.000000	40.250000
max	698.000000	3701.000000

	yhat	yhat_upper	yhat_lower
DATE
2024-03-23	93763.73	110608.42	75509.07
2024-03-30	98443.33	116404.55	81232.11
2024-04-06	108617.99	126353.05	92715.61
2024-04-13	115180.50	130721.01	97186.91

	ds	trend	yhat_lower	yhat_upper	trend_lower	trend_upper	additive_terms	additive_terms_lower	additive_terms_upper	yearly	yearly_lower	yearly_upper	yhat
0	2023-03-11	19485.394281	15920.064503	19464.224109	19483.843303	19488.275912	-1759.759261	-1759.759261	-1759.759261	-1759.759261	-1759.759261	-1759.759261	17725.635020
1	2023-03-18	19247.662322	16487.443642	19995.816392	19232.930247	19265.132350	-1028.066953	-1028.066953	-1028.066953	-1028.066953	-1028.066953	-1028.066953	18219.595369
2	2023-03-25	19009.930363	17453.491002	20805.587420	18976.524053	19051.408845	49.613071	49.613071	49.613071	49.613071	49.613071	49.613071	19059.543434
3	2023-04-01	18772.198405	18148.134072	21589.310780	18712.260886	18844.356146	1040.777482	1040.777482	1040.777482	1040.777482	1040.777482	1040.777482	19812.975887
4	2023-04-08	18534.466446	18501.014865	21900.768394	18440.762472	18635.148118	1679.434987	1679.434987	1679.434987	1679.434987	1679.434987	1679.434987	20213.901433

	Percentile	Value
0	75	40.25
1	80	52.00
2	85	72.00
3	90	115.00
4	95	231.00
5	96	299.00
6	97	402.00
7	98	594.90
8	99	1612.80
9	100	3701.00

	ds	yhat	yhat_lower	yhat_upper
195	2024-09-22	9112.684434	5202.013051	13014.200778
196	2024-09-29	8509.657160	4758.182902	12584.130980
197	2024-10-06	7670.426502	3965.991301	11509.779202
198	2024-10-13	6735.039638	2736.616109	10611.176965
199	2024-10-20	6006.095057	1857.985922	10155.324535

	ds	yhat
148	2023-10-29	13004.568376
149	2023-11-05	12468.837648
150	2023-11-12	11637.278663
151	2023-11-19	10663.885980
152	2023-11-26	9845.748470

	ds	trend	yhat_lower	yhat_upper	trend_lower	trend_upper	additive_terms	additive_terms_lower	additive_terms_upper	yearly	yearly_lower	yearly_upper	multiplicative_terms	multiplicative_terms_lower	multiplicative_terms_upper	yhat
135	2023-07-02	-110.053709	-156.706731	201.795569	-110.053709	-110.053709	135.114942	135.114942	135.114942	135.114942	135.114942	135.114942	0.0	0.0	0.0	25.061234
136	2023-07-09	-114.546123	-197.287970	165.726691	-114.818539	-114.317650	100.385351	100.385351	100.385351	100.385351	100.385351	100.385351	0.0	0.0	0.0	-14.160772

	ds	y
143	2023-09-30	4307.0
144	2023-10-07	4330.0
145	2023-10-14	4467.0
146	2023-10-21	4254.0
147	2023-10-28	3974.0

	count	mean	std	min	25%	50%	75%	max
State
CO	4409784.0	158.650056	426.420239	0.04	11.0	43.0	133.0	11219.0
OR	2645953.0	156.398806	356.790827	0.04	13.0	49.0	144.0	11100.0
UT	1460597.0	196.563396	471.057990	0.04	13.0	51.0	167.0	14773.0
WA	3049667.0	145.675230	331.236478	0.04	12.0	48.0	138.0	9243.0

	Region	Filter	Best 13 Weeks (Forecasted)	Absolute Value (Unit Sales)	Absolute Value Percentage for total weeks
0	NorthWest	Brand + Category	29 - 41	300731	30.62%
1	Central	Brand + Category	25 - 37	272066	31.24%
2	NorthWest	Flavor + Caloric + Category	30 - 42	173096	30.98%
3	Central	Flavor + Caloric + Category	32 - 44	108107	31.14%

	yhat	yhat_upper	yhat_lower
DATE
2023-09-02	285476.40	298778.29	271838.69
2023-09-09	288354.50	300856.25	276032.32
2023-09-16	289057.31	302472.57	274802.73
2023-09-23	287248.53	300958.92	273797.52
2023-09-30	285106.30	298506.32	271294.70
2023-10-07	284271.78	298363.42	271306.87
2023-10-14	283922.05	297361.37	270622.10
2023-10-21	282352.50	296781.29	268788.65
2023-10-28	279677.03	294189.70	266454.82
2023-11-04	277865.93	290756.87	264269.66
2023-11-11	278087.04	290981.57	264992.78
2023-11-18	279119.81	293013.24	265333.38
2023-11-25	279277.11	293360.46	265840.58

	DATE	UNIT_SALES
0	2020-12-05	1961.0
1	2020-12-12	1625.0

	Year	Unit_Sales
0	2020	4187.0
1	2021	50501.0
2	2022	7165.0
3	2023	156.0

	Year	Unit_Sales
0	2023	25.061234
1	2024	0.000000

	Year	Unit_Sales
0	2023	44463.494491
1	2024	173997.204651

Swire-CC Innovation Product Forecasting - Modelling

Authors: ¶

Table of Contents

1. Introduction¶

2. Project Goal & Analytics Approach

3. Data Description & Missing Value Analysis

4. Q1 Modelling¶

5. Q2 Modelling¶

6. Q3 Modelling¶

7. Q6 Modelling¶

8. Q4 Modelling¶

9. Q5 Modelling¶

10. Q7 Modelling¶

11. Team Contribution¶

Business Problem Introduction ¶

Innovative Product Business ¶

Business Problem ¶

Project Goal ¶

Modelling Notebook Purpose¶

Modelling Approach & Assumptions

Approach

Filtering Strategy

Assumptions

Data Preparation¶

Forecasting Models Considered¶

Q1

Prophet Model

Data Preparation for Modelling

Model Validation¶

Error Metrics¶

PLUM + SWIRE + DIET

Accounting for remaining parameters (Category & Brand) - Historic Ratio Method

Summary

Q2 ¶

Data Preparation¶

Model Performance Metrics Function¶

Modeling Process

Arima¶

SARIMA¶

Prophet¶

Exponential Smoothing¶

Performance Summary¶

Predicting Easter 2024 using Prophet¶

Summary¶

Q3

Data Preparation¶

Kiwano + Diet + Swire¶

Prophet Model¶

Train and Test Split:¶

Validation¶

Demand Forecasting for upcoming one year:¶

Summary¶

Q4 ¶

Mulberry + Swire + Diet

Prophet Model¶

Swire + Mulberry + Diet¶

Only Mulberry Prophet model¶

Mulberry + Diet¶

Mulberry + Swire¶

Diet + Mulberry + Sparkling¶

Prophet Model for validation of Sparkling Water + Diet + Mulberry¶

Training on entire dataset for forecasting sales for the next one year¶

Q4 Summary¶

Q6 ¶

Data Preparation¶

Model Performance Metrics Function¶

ARIMA¶

SARIMA¶

Prophet¶

Exponential Smoothing¶

Performance Summary¶

Predicting 6 month using Exponential Smoothing¶

Summary¶

Model Performance Summary¶

Findings and Conclusion¶

Q5

Central vs NW Statistical Analysis:

Greetingle + ENH Water Filtered Modelling:

Woodsy Yellow + Regular + Enh Water Filter Modelling¶

Q7

1. Introduction ¶

4. Q1 Modelling ¶

5. Q2 Modelling ¶

6. Q3 Modelling ¶

7. Q6 Modelling ¶

8. Q4 Modelling ¶

9. Q5 Modelling ¶

10. Q7 Modelling ¶

11. Team Contribution ¶